Fork me on GitHub

Changes in / [3e4e196:eee976c2] in git


Ignore:
Files:
3 deleted
11 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    r3e4e196 reee976c2  
    635635tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf): \
    636636        external/Hector/H_VerticalQuadrupole.$(SrcSuf)
    637 tmp/external/TrackCovariance/AcceptanceClx.$(ObjSuf): \
    638         external/TrackCovariance/AcceptanceClx.$(SrcSuf)
    639637tmp/external/TrackCovariance/ObsTrk.$(ObjSuf): \
    640638        external/TrackCovariance/ObsTrk.$(SrcSuf)
    … …  
    11471145        tmp/external/Hector/H_VerticalKicker.$(ObjSuf) \
    11481146        tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf) \
    1149         tmp/external/TrackCovariance/AcceptanceClx.$(ObjSuf) \
    11501147        tmp/external/TrackCovariance/ObsTrk.$(ObjSuf) \
    11511148        tmp/external/TrackCovariance/SolGeom.$(ObjSuf) \

  • cards/delphes_card_IDEA.tcl

    r3e4e196 reee976c2  
    11####################################################################                                l
    2 # FCC-ee IDEA detector model
    3 #
    4 # Authors: Elisa Fontanesi, Lorenzo Pezzotti, Massimiliano Antonello
     2# FCC-ee IDEA detector model                                                                                     
     3#                                                                                                   
     4# Authors: Elisa Fontanesi, Lorenzo Pezzotti, Massimiliano Antonello                                 
    55# email: efontane@bo.infn.it,
    6 #        lorenzo.pezzotti01@universitadipavia.it,
    7 #        m.antonello@uninsubria.it,
    8 #####################################################################
    9 #
     6#        lorenzo.pezzotti01@universitadipavia.it,                                                   
     7#        m.antonello@uninsubria.it,                                                                 
     8#####################################################################                               
     9#       
    1010#######################################
    1111# Order of execution of various modules
    … …  
    1919  MuonTrackingEfficiency
    2020
    21   TrackMergerPre
    22   TrackSmearing
    23 
    24   TrackMerger
     21  ChargedHadronMomentumSmearing
     22  ElectronMomentumSmearing
     23  MuonMomentumSmearing
     24
     25  TrackMerger
    2526  Calorimeter
    2627  EFlowMerger
    … …  
    3334  ElectronIsolation
    3435
    35   MuonFilter
    3636  MuonEfficiency
    3737  MuonIsolation
    … …  
    4242  GenJetFinder
    4343  GenMissingET
    44 
     44 
    4545  FastJetFinder
    4646
    … …  
    9898}
    9999
    100 #       (pt <= 0.1)                                     * (0.00) +
     100#       (pt <= 0.1)                                     * (0.00) + 
    101101#       (abs(eta) <= 3.0)               * (pt > 0.1)    * (1.00) +
    102102#       (abs(eta) > 3)                                  * (0.00)
    … …  
    119119        (energy < 0.5 && energy >= 0.3) * (abs(eta) <= 3.0)    * (0.65) +
    120120        (energy < 0.3) * (abs(eta) <= 3.0)                     * (0.06)
    121     }
     121    } 
    122122}
    123123
    … …  
    140140}
    141141
     142
     143########################################
     144# Momentum resolution for charged tracks
     145########################################
     146
     147module MomentumSmearing ChargedHadronMomentumSmearing {
     148    set InputArray ChargedHadronTrackingEfficiency/chargedHadrons
     149    set OutputArray chargedHadrons
     150
     151
     152    # Resolution given in dpT/pT.
     153    # IDEAdet
     154    set ResolutionFormula {
     155        (abs(eta) <= 3.0)                   * sqrt(0.0001145^2 + 0.0002024^2*pt + (pt*2.093e-005)^2)
     156    }
     157}
     158
     159###################################
     160# Momentum resolution for electrons
     161###################################
     162
     163module MomentumSmearing ElectronMomentumSmearing {
     164    set InputArray ElectronTrackingEfficiency/electrons
     165    set OutputArray electrons
     166
     167    # Resolution given in dpT/pT.
     168    # IDEAdet
     169    set ResolutionFormula {
     170        (abs(eta) <= 3.0)                   * sqrt(0.0001145^2 + 0.0002024^2*pt + (pt*2.093e-005)^2)
     171   } 
     172}
     173
     174###############################
     175# Momentum resolution for muons
     176###############################
     177
     178module MomentumSmearing MuonMomentumSmearing {
     179    set InputArray MuonTrackingEfficiency/muons
     180    set OutputArray muons
     181
     182    # Resolution given in dpT/pT.
     183    # IDEAdet
     184    set ResolutionFormula {
     185        (abs(eta) <= 3.0)                   * sqrt(0.0001145^2 + 0.0002024^2*pt + (pt*2.093e-005)^2)
     186    }
     187}
     188
    142189##############
    143190# Track merger
    144191##############
    145192
    146 module Merger TrackMergerPre {
    147 # add InputArray InputArray
    148   add InputArray ChargedHadronTrackingEfficiency/chargedHadrons
    149   add InputArray ElectronTrackingEfficiency/electrons
    150   add InputArray MuonTrackingEfficiency/muons
    151   set OutputArray tracks
    152 }
    153 
    154 
    155 ########################################
    156 # Smearing for charged tracks
    157 ########################################
    158 
    159 module TrackCovariance TrackSmearing {
    160     set InputArray TrackMergerPre/tracks
    161     set OutputArray tracks
    162 
    163     set Bz 2.0
    164 
    165     ## minimum number of hits to accept a track
    166     set NMinHits 6
    167 
    168     ## uses https://raw.githubusercontent.com/selvaggi/FastTrackCovariance/master/GeoIDEA_BASE.txt
    169     set DetectorGeometry {
    170 
    171 
    172       # Layer type 1 = R (barrel) or 2 = z (forward/backward)
    173       # Layer label
    174       # Minimum dimension z for barrel or R for forward
    175       # Maximum dimension z for barrel or R for forward
    176       # R/z location of layer
    177       # Thickness (meters)
    178       # Radiation length (meters)
    179       # Number of measurements in layers (1D or 2D)
    180       # Stereo angle (rad) - 0(pi/2) = axial(z) layer - Upper side
    181       # Stereo angle (rad) - 0(pi/2) = axial(z) layer - Lower side
    182       # Resolution Upper side (meters) - 0 = no measurement
    183       # Resolution Lower side (meters) - 0 = no measurement
    184       # measurement flag = T, scattering only = F
    185 
    186 
    187       # barrel  name       zmin   zmax   r        w (m)      X0        n_meas  th_up (rad) th_down (rad)    reso_up (m)   reso_down (m)  flag
    188 
    189       1        PIPE       -100    100    0.015    0.0012    0.35276    0        0          0                0             0              0
    190       1        VTXLOW     -0.12   0.12   0.017    0.00028   0.0937     2        0          1.5708           3e-006        3e-006         1
    191       1        VTXLOW     -0.16   0.16   0.023    0.00028   0.0937     2        0          1.5708           3e-006        3e-006         1
    192       1        VTXLOW     -0.16   0.16   0.031    0.00028   0.0937     2        0          1.5708           3e-006        3e-006         1
    193       1        VTXHIGH    -1      1      0.32     0.00047   0.0937     2        0          1.5708           7e-006        7e-006         1
    194       1        VTXHIGH    -1.05   1.05   0.34     0.00047   0.0937     2        0          1.5708           7e-006        7e-006         1
    195 
    196       # endcap  name       rmin   rmax   z        w (m)      X0        n_meas   th_up (rad)  th_down (rad)   reso_up (m)   reso_down (m) flag
    197 
    198       2        VTXDSK      0.141  0.3   -0.92     0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
    199       2        VTXDSK      0.138  0.3   -0.9      0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
    200       2        VTXDSK      0.065  0.3   -0.42     0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
    201       2        VTXDSK      0.062  0.3   -0.4      0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
    202       2        VTXDSK      0.062  0.3    0.4      0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
    203       2        VTXDSK      0.065  0.3    0.42     0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
    204       2        VTXDSK      0.138  0.3    0.9      0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
    205       2        VTXDSK      0.141  0.3    0.92     0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
    206 
    207       1 DCHCANI -2.125 2.125 0.345 0.0002 0.237223 0 0 0 0 0 0
    208       1 DCH -2 2 0.36 0.0147748 1400 1 0.0203738 0 0.0001 0 1
    209       1 DCH -2 2 0.374775 0.0147748 1400 1 -0.0212097 0 0.0001 0 1
    210       1 DCH -2 2 0.38955 0.0147748 1400 1 0.0220456 0 0.0001 0 1
    211       1 DCH -2 2 0.404324 0.0147748 1400 1 -0.0228814 0 0.0001 0 1
    212       1 DCH -2 2 0.419099 0.0147748 1400 1 0.0237172 0 0.0001 0 1
    213       1 DCH -2 2 0.433874 0.0147748 1400 1 -0.024553 0 0.0001 0 1
    214       1 DCH -2 2 0.448649 0.0147748 1400 1 0.0253888 0 0.0001 0 1
    215       1 DCH -2 2 0.463423 0.0147748 1400 1 -0.0262245 0 0.0001 0 1
    216       1 DCH -2 2 0.478198 0.0147748 1400 1 0.0270602 0 0.0001 0 1
    217       1 DCH -2 2 0.492973 0.0147748 1400 1 -0.0278958 0 0.0001 0 1
    218       1 DCH -2 2 0.507748 0.0147748 1400 1 0.0287314 0 0.0001 0 1
    219       1 DCH -2 2 0.522523 0.0147748 1400 1 -0.029567 0 0.0001 0 1
    220       1 DCH -2 2 0.537297 0.0147748 1400 1 0.0304025 0 0.0001 0 1
    221       1 DCH -2 2 0.552072 0.0147748 1400 1 -0.031238 0 0.0001 0 1
    222       1 DCH -2 2 0.566847 0.0147748 1400 1 0.0320734 0 0.0001 0 1
    223       1 DCH -2 2 0.581622 0.0147748 1400 1 -0.0329088 0 0.0001 0 1
    224       1 DCH -2 2 0.596396 0.0147748 1400 1 0.0337442 0 0.0001 0 1
    225       1 DCH -2 2 0.611171 0.0147748 1400 1 -0.0345795 0 0.0001 0 1
    226       1 DCH -2 2 0.625946 0.0147748 1400 1 0.0354147 0 0.0001 0 1
    227       1 DCH -2 2 0.640721 0.0147748 1400 1 -0.0362499 0 0.0001 0 1
    228       1 DCH -2 2 0.655495 0.0147748 1400 1 0.0370851 0 0.0001 0 1
    229       1 DCH -2 2 0.67027 0.0147748 1400 1 -0.0379202 0 0.0001 0 1
    230       1 DCH -2 2 0.685045 0.0147748 1400 1 0.0387552 0 0.0001 0 1
    231       1 DCH -2 2 0.69982 0.0147748 1400 1 -0.0395902 0 0.0001 0 1
    232       1 DCH -2 2 0.714595 0.0147748 1400 1 0.0404252 0 0.0001 0 1
    233       1 DCH -2 2 0.729369 0.0147748 1400 1 -0.04126 0 0.0001 0 1
    234       1 DCH -2 2 0.744144 0.0147748 1400 1 0.0420949 0 0.0001 0 1
    235       1 DCH -2 2 0.758919 0.0147748 1400 1 -0.0429296 0 0.0001 0 1
    236       1 DCH -2 2 0.773694 0.0147748 1400 1 0.0437643 0 0.0001 0 1
    237       1 DCH -2 2 0.788468 0.0147748 1400 1 -0.044599 0 0.0001 0 1
    238       1 DCH -2 2 0.803243 0.0147748 1400 1 0.0454336 0 0.0001 0 1
    239       1 DCH -2 2 0.818018 0.0147748 1400 1 -0.0462681 0 0.0001 0 1
    240       1 DCH -2 2 0.832793 0.0147748 1400 1 0.0471025 0 0.0001 0 1
    241       1 DCH -2 2 0.847568 0.0147748 1400 1 -0.0479369 0 0.0001 0 1
    242       1 DCH -2 2 0.862342 0.0147748 1400 1 0.0487713 0 0.0001 0 1
    243       1 DCH -2 2 0.877117 0.0147748 1400 1 -0.0496055 0 0.0001 0 1
    244       1 DCH -2 2 0.891892 0.0147748 1400 1 0.0504397 0 0.0001 0 1
    245       1 DCH -2 2 0.906667 0.0147748 1400 1 -0.0512738 0 0.0001 0 1
    246       1 DCH -2 2 0.921441 0.0147748 1400 1 0.0521079 0 0.0001 0 1
    247       1 DCH -2 2 0.936216 0.0147748 1400 1 -0.0529418 0 0.0001 0 1
    248       1 DCH -2 2 0.950991 0.0147748 1400 1 0.0537757 0 0.0001 0 1
    249       1 DCH -2 2 0.965766 0.0147748 1400 1 -0.0546095 0 0.0001 0 1
    250       1 DCH -2 2 0.980541 0.0147748 1400 1 0.0554433 0 0.0001 0 1
    251       1 DCH -2 2 0.995315 0.0147748 1400 1 -0.056277 0 0.0001 0 1
    252       1 DCH -2 2 1.01009 0.0147748 1400 1 0.0571106 0 0.0001 0 1
    253       1 DCH -2 2 1.02486 0.0147748 1400 1 -0.0579441 0 0.0001 0 1
    254       1 DCH -2 2 1.03964 0.0147748 1400 1 0.0587775 0 0.0001 0 1
    255       1 DCH -2 2 1.05441 0.0147748 1400 1 -0.0596108 0 0.0001 0 1
    256       1 DCH -2 2 1.06919 0.0147748 1400 1 0.0604441 0 0.0001 0 1
    257       1 DCH -2 2 1.08396 0.0147748 1400 1 -0.0612773 0 0.0001 0 1
    258       1 DCH -2 2 1.09874 0.0147748 1400 1 0.0621104 0 0.0001 0 1
    259       1 DCH -2 2 1.11351 0.0147748 1400 1 -0.0629434 0 0.0001 0 1
    260       1 DCH -2 2 1.12829 0.0147748 1400 1 0.0637763 0 0.0001 0 1
    261       1 DCH -2 2 1.14306 0.0147748 1400 1 -0.0646092 0 0.0001 0 1
    262       1 DCH -2 2 1.15784 0.0147748 1400 1 0.0654419 0 0.0001 0 1
    263       1 DCH -2 2 1.17261 0.0147748 1400 1 -0.0662746 0 0.0001 0 1
    264       1 DCH -2 2 1.18739 0.0147748 1400 1 0.0671071 0 0.0001 0 1
    265       1 DCH -2 2 1.20216 0.0147748 1400 1 -0.0679396 0 0.0001 0 1
    266       1 DCH -2 2 1.21694 0.0147748 1400 1 0.068772 0 0.0001 0 1
    267       1 DCH -2 2 1.23171 0.0147748 1400 1 -0.0696042 0 0.0001 0 1
    268       1 DCH -2 2 1.24649 0.0147748 1400 1 0.0704364 0 0.0001 0 1
    269       1 DCH -2 2 1.26126 0.0147748 1400 1 -0.0712685 0 0.0001 0 1
    270       1 DCH -2 2 1.27604 0.0147748 1400 1 0.0721005 0 0.0001 0 1
    271       1 DCH -2 2 1.29081 0.0147748 1400 1 -0.0729324 0 0.0001 0 1
    272       1 DCH -2 2 1.30559 0.0147748 1400 1 0.0737642 0 0.0001 0 1
    273       1 DCH -2 2 1.32036 0.0147748 1400 1 -0.0745958 0 0.0001 0 1
    274       1 DCH -2 2 1.33514 0.0147748 1400 1 0.0754274 0 0.0001 0 1
    275       1 DCH -2 2 1.34991 0.0147748 1400 1 -0.0762589 0 0.0001 0 1
    276       1 DCH -2 2 1.36468 0.0147748 1400 1 0.0770903 0 0.0001 0 1
    277       1 DCH -2 2 1.37946 0.0147748 1400 1 -0.0779215 0 0.0001 0 1
    278       1 DCH -2 2 1.39423 0.0147748 1400 1 0.0787527 0 0.0001 0 1
    279       1 DCH -2 2 1.40901 0.0147748 1400 1 -0.0795837 0 0.0001 0 1
    280       1 DCH -2 2 1.42378 0.0147748 1400 1 0.0804147 0 0.0001 0 1
    281       1 DCH -2 2 1.43856 0.0147748 1400 1 -0.0812455 0 0.0001 0 1
    282       1 DCH -2 2 1.45333 0.0147748 1400 1 0.0820762 0 0.0001 0 1
    283       1 DCH -2 2 1.46811 0.0147748 1400 1 -0.0829068 0 0.0001 0 1
    284       1 DCH -2 2 1.48288 0.0147748 1400 1 0.0837373 0 0.0001 0 1
    285       1 DCH -2 2 1.49766 0.0147748 1400 1 -0.0845677 0 0.0001 0 1
    286       1 DCH -2 2 1.51243 0.0147748 1400 1 0.0853979 0 0.0001 0 1
    287       1 DCH -2 2 1.52721 0.0147748 1400 1 -0.086228 0 0.0001 0 1
    288       1 DCH -2 2 1.54198 0.0147748 1400 1 0.087058 0 0.0001 0 1
    289       1 DCH -2 2 1.55676 0.0147748 1400 1 -0.0878879 0 0.0001 0 1
    290       1 DCH -2 2 1.57153 0.0147748 1400 1 0.0887177 0 0.0001 0 1
    291       1 DCH -2 2 1.58631 0.0147748 1400 1 -0.0895474 0 0.0001 0 1
    292       1 DCH -2 2 1.60108 0.0147748 1400 1 0.0903769 0 0.0001 0 1
    293       1 DCH -2 2 1.61586 0.0147748 1400 1 -0.0912063 0 0.0001 0 1
    294       1 DCH -2 2 1.63063 0.0147748 1400 1 0.0920356 0 0.0001 0 1
    295       1 DCH -2 2 1.64541 0.0147748 1400 1 -0.0928647 0 0.0001 0 1
    296       1 DCH -2 2 1.66018 0.0147748 1400 1 0.0936937 0 0.0001 0 1
    297       1 DCH -2 2 1.67495 0.0147748 1400 1 -0.0945226 0 0.0001 0 1
    298       1 DCH -2 2 1.68973 0.0147748 1400 1 0.0953514 0 0.0001 0 1
    299       1 DCH -2 2 1.7045 0.0147748 1400 1 -0.09618 0 0.0001 0 1
    300       1 DCH -2 2 1.71928 0.0147748 1400 1 0.0970085 0 0.0001 0 1
    301       1 DCH -2 2 1.73405 0.0147748 1400 1 -0.0978369 0 0.0001 0 1
    302       1 DCH -2 2 1.74883 0.0147748 1400 1 0.0986651 0 0.0001 0 1
    303       1 DCH -2 2 1.7636 0.0147748 1400 1 -0.0994932 0 0.0001 0 1
    304       1 DCH -2 2 1.77838 0.0147748 1400 1 0.100321 0 0.0001 0 1
    305       1 DCH -2 2 1.79315 0.0147748 1400 1 -0.101149 0 0.0001 0 1
    306       1 DCH -2 2 1.80793 0.0147748 1400 1 0.101977 0 0.0001 0 1
    307       1 DCH -2 2 1.8227 0.0147748 1400 1 -0.102804 0 0.0001 0 1
    308       1 DCH -2 2 1.83748 0.0147748 1400 1 0.103632 0 0.0001 0 1
    309       1 DCH -2 2 1.85225 0.0147748 1400 1 -0.104459 0 0.0001 0 1
    310       1 DCH -2 2 1.86703 0.0147748 1400 1 0.105286 0 0.0001 0 1
    311       1 DCH -2 2 1.8818 0.0147748 1400 1 -0.106113 0 0.0001 0 1
    312       1 DCH -2 2 1.89658 0.0147748 1400 1 0.10694 0 0.0001 0 1
    313       1 DCH -2 2 1.91135 0.0147748 1400 1 -0.107766 0 0.0001 0 1
    314       1 DCH -2 2 1.92613 0.0147748 1400 1 0.108593 0 0.0001 0 1
    315       1 DCH -2 2 1.9409 0.0147748 1400 1 -0.109419 0 0.0001 0 1
    316       1 DCH -2 2 1.95568 0.0147748 1400 1 0.110246 0 0.0001 0 1
    317       1 DCH -2 2 1.97045 0.0147748 1400 1 -0.111072 0 0.0001 0 1
    318       1 DCH -2 2 1.98523 0.0147748 1400 1 0.111898 0 0.0001 0 1
    319       1 DCH -2 2 2 0.0147748 1400 1 -0.112723 0 0.0001 0 1
    320       1 DCHCANO -2.125 2.125 2.02 0.02 1.667 0 0 0 0 0 0
    321       1 BSILWRP -2.35 2.35 2.04 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    322       1 BSILWRP -2.35 2.35 2.06 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    323       1 MAG -2.5 2.5 2.25 0.05 0.0658 0 0 0 0 0 0
    324       1 BPRESH -2.55 2.55 2.45 0.02 1 2 0 1.5708 7e-005 0.01 1
    325 
    326 
    327       2 DCHWALL 0.345 2.02 2.125 0.25 5.55 0 0 0 0 0 0
    328       2 DCHWALL 0.345 2.02 -2.125 0.25 5.55 0 0 0 0 0 0
    329       2 FSILWRP 0.354 2.02 -2.32 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    330       2 FSILWRP 0.35 2.02 -2.3 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    331       2 FSILWRP 0.35 2.02 2.3 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    332       2 FSILWRP 0.354 2.02 2.32 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    333       2 FRAD 0.38 2.09 2.49 0.0043 0.005612 0 0 0 0 0 0
    334       2 FRAD 0.38 2.09 -2.49 0.0043 0.005612 0 0 0 0 0 0
    335       2 FPRESH 0.39 2.43 -2.55 0.02 1 2 0 1.5708 7e-005 0.01 1
    336       2 FPRESH 0.39 2.43 2.55 0.02 1 2 0 1.5708 7e-005 0.01 1
    337     }
    338 
    339 }
    340 
    341 
    342 ##############
    343 # Track merger
    344 ##############
    345 
    346193module Merger TrackMerger {
    347194# add InputArray InputArray
    348   add InputArray TrackSmearing/tracks
     195  add InputArray ChargedHadronMomentumSmearing/chargedHadrons
     196  add InputArray ElectronMomentumSmearing/electrons
     197  add InputArray MuonMomentumSmearing/muons
    349198  set OutputArray tracks
    350199}
    351200
    352201
    353 #############
    354 # Calorimeter
    355 #############
     202#############                                                                                                                         
     203# Calorimeter                                                                                                                                           
     204#############                                                                                                                                           
    356205module DualReadoutCalorimeter Calorimeter {
    357206  set ParticleInputArray ParticlePropagator/stableParticles
    … …  
    375224    set pi [expr {acos(-1)}]
    376225
    377     # Lists of the edges of each tower in eta and phi;
    378     # each list starts with the lower edge of the first tower;
    379     # the list ends with the higher edged of the last tower.
    380     # Barrel:  deta=0.02 towers up to |eta| <= 0.88 ( up to 45°)
    381     # Endcaps: deta=0.02 towers up to |eta| <= 3.0 (8.6° = 100 mrad)
    382     # Cell size: about 6 cm x 6 cm
    383 
    384     #barrel:
     226    # Lists of the edges of each tower in eta and phi;                                                                                           
     227    # each list starts with the lower edge of the first tower;                                                                                                   
     228    # the list ends with the higher edged of the last tower.                                                                                       
     229    # Barrel:  deta=0.02 towers up to |eta| <= 0.88 ( up to 45°)                                                                                       
     230    # Endcaps: deta=0.02 towers up to |eta| <= 3.0 (8.6° = 100 mrad)                                                                                           
     231    # Cell size: about 6 cm x 6 cm                                                 
     232
     233    #barrel:                                                                                       
    385234    set PhiBins {}
    386235    for {set i -120} {$i <= 120} {incr i} {
    387236        add PhiBins [expr {$i * $pi/120}]
    388237    }
    389     #deta=0.02 units for |eta| <= 0.88
     238    #deta=0.02 units for |eta| <= 0.88                                                             
    390239    for {set i -44} {$i < 45} {incr i} {
    391240        set eta [expr {$i * 0.02}]
    … …  
    393242    }
    394243
    395     #endcaps:
     244    #endcaps:                                                                                       
    396245    set PhiBins {}
    397246    for {set i -120} {$i <= 120} {incr i} {
    398247        add PhiBins [expr {$i* $pi/120}]
    399248    }
    400     #deta=0.02 units for 0.88 < |eta| <= 3.0
    401     #first, from -3.0 to -0.88
     249    #deta=0.02 units for 0.88 < |eta| <= 3.0                                                       
     250    #first, from -3.0 to -0.88                                                                     
    402251    for {set i 1} {$i <=106} {incr i} {
    403252        set eta [expr {-3.00 + $i * 0.02}]
    404253        add EtaPhiBins $eta $PhiBins
    405254    }
    406     #same for 0.88 to 3.0
     255    #same for 0.88 to 3.0                                                                           
    407256    for  {set i 1} {$i <=106} {incr i} {
    408257        set eta [expr {0.88 + $i * 0.02}]
    … …  
    410259    }
    411260
    412     # default energy fractions {abs(PDG code)} {Fecal Fhcal}
     261    # default energy fractions {abs(PDG code)} {Fecal Fhcal}                                                                                 
    413262    add EnergyFraction {0} {0.0 1.0}
    414     # energy fractions for e, gamma and pi0
     263    # energy fractions for e, gamma and pi0                                                                                                           
    415264    add EnergyFraction {11} {1.0 0.0}
    416265    add EnergyFraction {22} {1.0 0.0}
    417266    add EnergyFraction {111} {1.0 0.0}
    418     # energy fractions for muon, neutrinos and neutralinos
     267    # energy fractions for muon, neutrinos and neutralinos                                                                             
    419268    add EnergyFraction {12} {0.0 0.0}
    420269    add EnergyFraction {13} {0.0 0.0}
    … …  
    426275    add EnergyFraction {1000035} {0.0 0.0}
    427276    add EnergyFraction {1000045} {0.0 0.0}
    428     # energy fractions for K0short and Lambda
     277    # energy fractions for K0short and Lambda                                                                                                           
    429278    add EnergyFraction {310} {0.3 0.7}
    430279    add EnergyFraction {3122} {0.3 0.7}
    431280
    432281
    433     # set ECalResolutionFormula {resolution formula as a function of eta and energy}
     282    # set ECalResolutionFormula {resolution formula as a function of eta and energy}                               
    434283    set ECalResolutionFormula {
    435284    (abs(eta) <= 0.88 )                     * sqrt(energy^2*0.01^2 + energy*0.11^2)+
    … …  
    437286    }
    438287
    439     # set HCalResolutionFormula {resolution formula as a function of eta and energy}
     288    # set HCalResolutionFormula {resolution formula as a function of eta and energy}                                               
    440289    set HCalResolutionFormula {
    441290    (abs(eta) <= 0.88 )                     * sqrt(energy^2*0.01^2 + energy*0.30^2)+
    … …  
    503352}
    504353
    505 
    506 #################
    507 # Muon filter
    508 #################
    509 
    510 module PdgCodeFilter MuonFilter {
    511   set InputArray Calorimeter/eflowTracks
    512   set OutputArray muons
    513   set Invert true
    514   add PdgCode {13}
    515   add PdgCode {-13}
    516 }
    517 
    518 
    519354#####################
    520355# Electron efficiency
    … …  
    528363
    529364  # efficiency formula for electrons
    530   set EfficiencyFormula {
     365  set EfficiencyFormula {         
    531366        (energy < 2.0)                                         * (0.000)+
    532367        (energy >= 2.0) * (abs(eta) <= 0.88)                   * (0.99) +
    … …  
    553388}
    554389
    555 
    556390#################
    557391# Muon efficiency
    … …  
    559393
    560394module Efficiency MuonEfficiency {
    561   set InputArray MuonFilter/muons
     395  set InputArray MuonMomentumSmearing/muons
    562396  set OutputArray muons
    563397
    … …  
    565399
    566400  # efficiency formula for muons
    567   set EfficiencyFormula {
     401  set EfficiencyFormula {                                   
    568402        (energy < 2.0)                                         * (0.000)+
    569403        (energy >= 2.0) * (abs(eta) <= 0.88)                   * (0.99) +
    … …  
    713547
    714548  # add EfficiencyFormula {abs(PDG code)} {efficiency formula as a function of eta and pt}
    715 
     549 
    716550  # default efficiency formula (misidentification rate)
    717551  add EfficiencyFormula {0} {0.01}
    … …  
    769603module TreeWriter TreeWriter {
    770604    # add Branch InputArray BranchName BranchClass
    771 
     605   
    772606    add Branch Delphes/allParticles Particle GenParticle
    773607
    774608    add Branch TrackMerger/tracks Track Track
    775609    add Branch Calorimeter/towers Tower Tower
    776 
     610   
    777611    add Branch Calorimeter/eflowTracks EFlowTrack Track
    778612    add Branch Calorimeter/eflowPhotons EFlowPhoton Tower
    … …  
    782616    add Branch PhotonEfficiency/photons PhotonEff Photon
    783617    add Branch PhotonIsolation/photons PhotonIso Photon
    784 
     618   
    785619    add Branch GenJetFinder/jets GenJet Jet
    786620    add Branch GenMissingET/momentum GenMissingET MissingET
    787 
     621   
    788622    add Branch UniqueObjectFinder/jets Jet Jet
    789623    add Branch UniqueObjectFinder/electrons Electron Electron
    790624    add Branch UniqueObjectFinder/photons Photon Photon
    791625    add Branch UniqueObjectFinder/muons Muon Muon
    792 
    793     add Branch JetEnergyScale/jets AntiKtJet Jet
    794 
     626   
     627    add Branch JetEnergyScale/jets AntiKtJet Jet 
     628   
    795629    add Branch MissingET/momentum MissingET MissingET
    796630    add Branch ScalarHT/energy ScalarHT ScalarHT
    797631}
     632

  • external/TrackCovariance/ObsTrk.cc

    r3e4e196 reee976c2  
    4444        fCovILC = CovToILC(fCov);
    4545}
    46 
    47 // x[3] track origin, p[3] track momentum at origin, Q charge, B magnetic field in Tesla
    48 ObsTrk::ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC)
    49 {
    50         fGC = GC;
    51         fGenX.SetXYZ(x[0],x[1],x[2]);
    52         fGenP.SetXYZ(p[0],p[1],p[2]);
    53         fGenQ = Q;
    54         fB = B;
    55         fGenPar.ResizeTo(5);
    56         fGenParACTS.ResizeTo(6);
    57         fGenParILC.ResizeTo(5);
    58         fObsPar.ResizeTo(5);
    59         fObsParACTS.ResizeTo(6);
    60         fObsParILC.ResizeTo(5);
    61         fCov.ResizeTo(5, 5);
    62         fCovACTS.ResizeTo(6, 6);
    63         fCovILC.ResizeTo(5, 5);
    64         fGenPar = XPtoPar(fGenX, fGenP, Q);
    65         fGenParACTS = ParToACTS(fGenPar);
    66         fGenParILC = ParToILC(fGenPar);
    67         /*
    68         cout << "ObsTrk::ObsTrk: fGenPar";
    69         for (Int_t i = 0; i < 5; i++)cout << fGenPar(i) << ", ";
    70         cout << endl;
    71         */
    72         fObsPar = GenToObsPar(fGenPar, fGC);
    73         fObsParACTS = ParToACTS(fObsPar);
    74         fObsParILC = ParToILC(fObsPar);
    75         fObsX = ParToX(fObsPar);
    76         fObsP = ParToP(fObsPar);
    77         fObsQ = ParToQ(fObsPar);
    78         fCovACTS = CovToACTS(fCov);
    79         fCovILC = CovToILC(fCov);
    80 }
    81 
    82 
    8346//
    8447// Destructor
    … …  
    13598        //
    13699        TVector3 Xval;
    137         Xval(0) = -D*TMath::Sin(phi0);
    138         Xval(1) =  D*TMath::Cos(phi0);
     100        Xval(0) = -D*TMath::Sin(phi0); 
     101        Xval(1) =  D*TMath::Cos(phi0); 
    139102        Xval(2) =  z0;
    140103        //
    … …  
    172135        // Check ranges
    173136        Double_t minPt = GC->GetMinPt ();
    174         //if (pt < minPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << std::endl;
     137        if (pt < minPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << std::endl;
    175138        Double_t maxPt = GC->GetMaxPt();
    176         //if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl;
     139        if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl;
    177140        Double_t minAn = GC->GetMinAng();
    178   //if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
    179         //      << " is below grid range of " << minAn << std::endl;
     141        if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
     142                << " is below grid range of " << minAn << std::endl;
    180143        Double_t maxAn = GC->GetMaxAng();
    181         //if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
    182         //      << " is above grid range of " << maxAn << std::endl;
     144        if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
     145                << " is above grid range of " << maxAn << std::endl;
    183146        //
    184147        TMatrixDSym Cov = GC->GetCov(pt, angd);
    … …  
    216179        pACTS(1) = 1000 * Par(3);       // z0 from m to mm
    217180        pACTS(2) = Par(1);                      // Phi0 is unchanged
    218         pACTS(3) = TMath::ATan2(1.0,Par(4));            // Theta in [0, pi] range
     181        pACTS(3) = TMath::ATan(1.0 / Par(4)) + TMath::PiOver2();                // Theta in [0, pi] range
    219182        pACTS(4) = Par(2) / (b*TMath::Sqrt(1 + Par(4)*Par(4)));         // q/p in GeV
    220183        pACTS(5) = 0.0;                         // Time: currently undefined
    … …  
    284247        return cILC;
    285248}
     249
     250
     251
     252       

  • external/TrackCovariance/ObsTrk.h

    r3e4e196 reee976c2  
    1919        // Prefix Gen marks variables before resolution smearing
    2020        //
    21 private:
     21private:       
    2222        Double_t fB;                                            // Solenoid magnetic field
    2323        SolGridCov *fGC;                                        // Covariance matrix grid
    … …  
    2525        Double_t fObsQ;                                 // Observed  track charge
    2626        TVector3 fGenX;                                 // Generated track origin (x,y,z)
    27         TVector3 fObsX;                                 // Observed  track origin (x,y,z) @ track min. approach
     27        TVector3 fObsX;                                 // Observed  track origin (x,y,z) @ track min. approach 
    2828        TVector3 fGenP;                                 // Generated track momentum at track origin
    2929        TVector3 fObsP;                                 // Observed  track momentum @ track minimum approach
    … …  
    4343        //
    4444        TVectorD ParToACTS(TVectorD Par);               // Parameter conversion
    45         TMatrixDSym CovToACTS(TMatrixDSym Cov); // Covariance
     45        TMatrixDSym CovToACTS(TMatrixDSym Cov); // Covariance 
    4646        //
    4747        // Conversion to ILC parametrization
    … …  
    5555        // x(3) track origin, p(3) track momentum at origin, Q charge, B magnetic field in Tesla
    5656        ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC); // Initialize and generate smeared track
    57         ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC);       // Initialize and generate smeared track
    58   // Destructor
     57        // Destructor
    5958        ~ObsTrk();
    6059        //

  • external/TrackCovariance/SolGridCov.cc

    r3e4e196 reee976c2  
    33#include <TMath.h>
    44#include <TVectorD.h>
    5 #include <TVector3.h>
    65#include <TMatrixDSym.h>
    76#include <TDecompChol.h>
    … …  
    3736{
    3837  delete[] fCov;
    39   delete fAcc;
    4038}
    4139
    … …  
    6058    }
    6159  }
    62 
    63 // Now make acceptance
    64 fAcc = new AcceptanceClx(G);
    6560}
    66 
    67 
    68 //
    69 Bool_t SolGridCov::IsAccepted(Double_t pt, Double_t Theta)
    70 {
    71         //
    72         // pt in GeV, Theta in degrees
    73         //
    74         Bool_t Accept = kFALSE;
    75         if (fAcc->HitNumber(pt, Theta) >= fNminHits)Accept = kTRUE;
    76         //
    77         return Accept;
    78 }
    79 //
    80 Bool_t SolGridCov::IsAccepted(Double_t *p)
    81 {
    82         //
    83         // pt in GeV, Theta in degrees
    84         //
    85         Bool_t Accept = kFALSE;
    86         Double_t pt = TMath::Sqrt(p[0] * p[0] + p[1] * p[1]);
    87         Double_t th = 180. * TMath::ATan2(pt, p[2])/TMath::Pi();
    88         if (fAcc->HitNumber(pt,th) >= fNminHits)Accept = kTRUE;
    89         //
    90         return Accept;
    91 }
    92 //
    93 Bool_t SolGridCov::IsAccepted(TVector3 p)
    94 {
    95         //
    96         // pt in GeV, Theta in degrees
    97         //
    98         Bool_t Accept = kFALSE;
    99         Double_t pt = p.Pt();
    100         Double_t th = 180.*TMath::ACos(p.CosTheta())/TMath::Pi();
    101         if (fAcc->HitNumber(pt,th) >= fNminHits)Accept = kTRUE;
    102         //
    103         return Accept;
    104 }
    105 
    106 
    10761// Find bin in grid
    10862Int_t SolGridCov::GetMinIndex(Double_t xval, Int_t N, TVectorD x)

  • external/TrackCovariance/SolGridCov.h

    r3e4e196 reee976c2  
    44#include <TVectorD.h>
    55#include <TMatrixDSym.h>
    6 #include "AcceptanceClx.h"
    76
    87class SolGeom;
    … …  
    1817  TVectorD fAnga;    // Array of angle points in degrees
    1918  TMatrixDSym *fCov; // Pointers to grid of covariance matrices
    20         AcceptanceClx *fAcc;                            // Pointer to acceptance class
    21         Int_t fNminHits;                                        // Minimum number of hits to accept track
    2219  // Service routines
    2320  Int_t GetMinIndex(Double_t xval, Int_t N, TVectorD x); // Find bin
    … …  
    3532  Double_t GetMaxAng() { return fAnga(fNang - 1); }
    3633  TMatrixDSym GetCov(Double_t pt, Double_t ang);
    37 
    38         // Acceptance related methods
    39         AcceptanceClx* AccPnt() { return fAcc; };                       // Return Acceptance class pointer
    40         void SetMinHits(Int_t MinHits) { fNminHits = MinHits; };        // Set minimum number of hits to accept (default = 6)
    41         Int_t GetMinHits() { return fNminHits; };
    42         Bool_t IsAccepted(Double_t pt, Double_t Theta);         // From pt (GeV) and theta (degrees)
    43         Bool_t IsAccepted(Double_t *p);                                         // From momentum components (GeV)
    44         Bool_t IsAccepted(TVector3 p);                                          // As above in Vector3 format
    45 
    4634};
    4735

  • external/TrackCovariance/SolTrack.cc

    r3e4e196 reee976c2  
    4545  fCov.ResizeTo(5, 5);
    4646}
    47 SolTrack::SolTrack(TVector3 x, TVector3 p, SolGeom* G)
    48 {
    49         fG = G;
    50         // Store momentum
    51         fp[0] = p(0); fp[1] = p(1); fp[2] = p(2);
    52         Double_t px = p(0); Double_t py = p(1); Double_t pz = p(2);
    53         fx[0] = x(0); fx[1] = x(1); fx[2] = x(2);
    54         Double_t xx = x(0); Double_t yy = x(1); Double_t zz = x(2);
    55         // Store parameters
    56         Double_t pt = TMath::Sqrt(px * px + py * py);
    57         Double_t Charge = 1.0;                                          // Don't worry about charge for now
    58         Double_t a = -Charge * G->B() * 0.2998;                 // Normalized speed of light
    59         Double_t C = a / (2 * pt);                                      // pt in GeV, B in Tesla, C in meters
    60         Double_t r2 = xx * xx + yy * yy;
    61         Double_t cross = xx * py - yy * px;
    62         Double_t T = TMath::Sqrt(pt * pt - 2 * a * cross + a * a * r2);
    63         Double_t phi0 = TMath::ATan2((py - a * xx) / T, (px + a * yy) / T);
    64         Double_t D;
    65         if (pt < 10.0) D = (T - pt) / a;
    66         else D = (-2 * cross + a * r2) / (T + pt);
    67         Double_t B = C * TMath::Sqrt((r2 - D * D) / (1 + 2 * C * D));
    68         Double_t st = TMath::ASin(B) / C;
    69         Double_t ct = pz / pt;
    70         Double_t z0 = zz - ct * st;
    71         fpar[0] = D;
    72         fpar[1] = phi0;
    73         fpar[2] = C;
    74         fpar[3] = z0;
    75         fpar[4] = ct;
    76         //cout << "SolTrack:: C = " << C << ", fpar[2] = " << fpar[2] << endl;
    77         //
    78         // Init covariances
    79         //
    80         fCov.ResizeTo(5, 5);
    81 }
    8247
    8348SolTrack::SolTrack(Double_t D, Double_t phi0, Double_t C, Double_t z0, Double_t ct, SolGeom *G)
    … …  
    10469SolTrack::~SolTrack()
    10570{
    106         delete[] & fp;
    107         delete[] & fpar;
    10871  fCov.Clear();
    10972}
    … …  
    169132  return kh;
    170133}
    171 
    172 // # of measurement layers hit
    173 Int_t SolTrack::nmHit()
    174 {
    175         Int_t kmh = 0;
    176         Double_t R; Double_t phi; Double_t zz;
    177         for (Int_t i = 0; i < fG->Nl(); i++)
    178         if (HitLayer(i, R, phi, zz))if (fG->isMeasure(i))kmh++;
    179         //
    180         return kmh;
    181 }
    182 //
    183 
    184 
    185134// List of layers hit with intersections
    186135Int_t SolTrack::HitList(Int_t *&ihh, Double_t *&rhh, Double_t *&zhh)
    … …  
    212161  return kmh;
    213162}
    214 
    215 // List of XYZ measurements without any error
    216 Int_t SolTrack::HitListXYZ(Int_t *&ihh, Double_t *&Xh, Double_t *&Yh, Double_t *&Zh)
    217 {
    218         //
    219         // Return lists of hits associated to a track for all measurement layers.
    220         // Return value is the total number of measurement hits
    221         // kmh = total number of measurement layers hit for given type
    222         // ihh = pointer to layer number
    223         // Xh, Yh, Zh = X, Y, Z of hit - No measurement error - No multiple scattering
    224         //
    225         //
    226         Int_t kmh = 0;  // Number of measurement layers hit
    227         for (Int_t i = 0; i < fG->Nl(); i++)
    228         {
    229                 Double_t R; Double_t phi; Double_t zz;
    230                 if (HitLayer(i, R, phi, zz)) // Only barrel type layers
    231                 {
    232                         if (fG->isMeasure(i))
    233                         {
    234                                 ihh[kmh] = i;
    235                                 Xh[kmh] = R*cos(phi);
    236                                 Yh[kmh] = R*sin(phi);
    237                                 Zh[kmh] = zz;
    238                                 kmh++;
    239                         }
    240                 }
    241         }
    242         //
    243         return kmh;
    244 }
    245 //
    246163// Covariance matrix estimation
    247 //
    248164void SolTrack::CovCalc(Bool_t Res, Bool_t MS)
    249165{
    250         //
    251         //
    252         // Input flags:
    253         //                              Res = .TRUE. turn on resolution effects/Use standard resolutions
    254         //                                        .FALSE. set all resolutions to 0
    255         //                              MS  = .TRUE. include Multiple Scattering
    256         //
    257         // Assumptions:
    258         // 1. Measurement layers can do one or two measurements
    259         // 2. On disks at constant z:
    260         //              - Upper side measurement is phi
    261         //              - Lower side measurement is R
    262         //
    263         // Fill list of layers hit
    264         //
    265         Int_t ntry = 0;
    266         Int_t ntrymax = 0;
    267         Int_t Nhit = nHit();                                    // Total number of layers hit
    268         Double_t *zhh = new Double_t[Nhit];             // z of hit
    269         Double_t *rhh = new Double_t[Nhit];             // r of hit
    270         Double_t *dhh = new Double_t[Nhit];             // distance of hit from origin
    271         Int_t    *ihh = new Int_t[Nhit];                // true index of layer
    272         Int_t kmh;                                                              // Number of measurement layers hit
    273         //
    274         kmh = HitList(ihh, rhh, zhh);                   // hit layer list
    275         Int_t mTot = 0;                                                 // Total number of measurements
    276         for (Int_t i = 0; i < Nhit; i++)
    277         {
    278                 dhh[i] = TMath::Sqrt(rhh[i] * rhh[i] + zhh[i] * zhh[i]);
    279                 if (fG->isMeasure(ihh[i])) mTot += fG->lND(ihh[i]);     // Count number of measurements
    280         }
    281         //
    282         // Order hit list by increasing distance from origin
    283         //
    284         Int_t    *hord = new Int_t[Nhit];               // hit order by increasing distance from origin
    285         TMath::Sort(Nhit, dhh, hord, kFALSE);   // Order by increasing distance from origin
    286         Double_t *zh = new Double_t[Nhit];              // d-ordered z of hit
    287         Double_t *rh = new Double_t[Nhit];              // d-ordered r of hit
    288         Int_t    *ih = new Int_t[Nhit];                 // d-ordered true index of layer
    289         for (Int_t i = 0; i < Nhit; i++)
    290         {
    291                 Int_t il = hord[i];                                     // Hit layer numbering
    292                 zh[i] = zhh[il];
    293                 rh[i] = rhh[il];
    294                 ih[i] = ihh[il];
    295         }
    296         //
    297         // Store interdistances and multiple scattering angles
    298         //
    299         Double_t sn2t = 1.0 / (1 + ct()*ct());                  //sin^2 theta of track
    300         Double_t cs2t = 1.0 - sn2t;                                             //cos^2 theta
    301         Double_t snt = TMath::Sqrt(sn2t);                               // sin theta
    302         Double_t cst = TMath::Sqrt(cs2t);                               // cos theta
    303         Double_t px0 = pt() * TMath::Cos(phi0());               // Momentum at minimum approach
    304         Double_t py0 = pt() * TMath::Sin(phi0());
    305         Double_t pz0 = pt() * ct();
    306         //
    307         TMatrixDSym dik(Nhit);  dik.Zero();             // Distances between layers
    308         Double_t *thms = new Double_t[Nhit];    // Scattering angles/plane
    309         Double_t *cs = new Double_t[Nhit];              // Cosine of angle with layer normal
    310         for (Int_t ii = 0; ii < Nhit; ii++)             // Hit layer loop
    311         {
    312                 Int_t i = ih[ii];                                       // Get true layer number
    313                 Double_t B = C()*TMath::Sqrt((rh[ii] * rh[ii] - D()*D()) / (1 + 2 * C()*D()));
    314                 Double_t pxi = px0*(1-2*B*B)-2*py0*B*TMath::Sqrt(1-B*B);                // Momentum at scattering layer
    315                 Double_t pyi = py0*(1-2*B*B)+2*px0*B*TMath::Sqrt(1-B*B);
    316                 Double_t pzi = pz0;
    317                 Double_t ArgRp = (rh[ii]*C() + (1 + C() * D())*D() / rh[ii]) / (1 + 2 * C()*D());
    318                 Double_t phi = phi0() + TMath::ASin(ArgRp);
    319                 Double_t nx = TMath::Cos(phi);          // Barrel layer normal
    320                 Double_t ny = TMath::Sin(phi);
    321                 Double_t nz = 0.0;
    322                 if (fG->lTyp(i) == 2)                                                           // this is Z layer
    323                 {
    324                         nx = 0.0;
    325                         ny = 0.0;
    326                         nz = 1.0;
    327                 }
    328                 Double_t corr = (pxi*nx + pyi * ny + pzi * nz) / p();
    329                 cs[ii] = corr;
    330                 Double_t Rlf = fG->lTh(i) / (corr*fG->lX0(i));          // Rad. length fraction
    331                 thms[ii] = 0.0136*TMath::Sqrt(Rlf)*(1.0 + 0.038*TMath::Log(Rlf)) / p();         // MS angle
    332                 if (!MS)thms[ii] = 0;
    333                 //
    334                 for (Int_t kk = 0; kk < ii; kk++)       // Fill distances between layers
    335                 {
    336                         //dik(ii, kk) = TMath::Sqrt(pow(rh[ii] - rh[kk], 2) + pow(zh[ii] - zh[kk], 2));
    337                         Double_t Ci = C();
    338                         dik(ii, kk) = (TMath::ASin(Ci*rh[ii])-TMath::ASin(Ci*rh[kk]))/(Ci*snt);
    339                         dik(kk, ii) = dik(ii, kk);
    340                 }
    341                 //
    342                 // Store momentum components for resolution correction cosines
    343                 //
    344                 Double_t *pRi = new Double_t[Nhit];
    345                 pRi[ii] = TMath::Abs(pxi * TMath::Cos(phi) + pyi * TMath::Sin(phi)); // Radial component
    346                 Double_t *pPhi = new Double_t[Nhit];
    347                 pPhi[ii] = TMath::Abs(pxi * TMath::Sin(phi) - pyi * TMath::Cos(phi)); // Phi component
    348         }
    349         //
    350         // Fill measurement covariance
    351         //
    352         Int_t *mTl = new Int_t[mTot];           // Pointer from measurement number to true layer number
    353         TMatrixDSym Sm(mTot); Sm.Zero();        // Measurement covariance
    354         TMatrixD Rm(mTot, 5);                           // Derivative matrix
    355         Int_t im = 0;
    356         //
    357         // Fill derivatives and error matrix with MS
    358         //
    359         Double_t AngMax = 90.; Double_t AngMx = AngMax * TMath::Pi() / 180.;
    360         Double_t csMin = TMath::Cos(AngMx);     // Set maximum angle wrt normal
    361         //
    362         for (Int_t ii = 0; ii < Nhit; ii++)
    363         {
    364                 Int_t i = ih[ii];                                       // True layer number
    365                 Int_t ityp  = fG->lTyp(i);                      // Layer type Barrel or Z
    366                 Int_t nmeai = fG->lND(i);                       // # measurements in layer
    367                 if (fG->isMeasure(i) && nmeai >0 && cs[ii] > csMin)
    368                 {
    369                         Double_t Di    = D();                           // Get true track parameters
    370                         Double_t phi0i = phi0();
    371                         Double_t Ci    = C();
    372                         Double_t z0i   = z0();
    373                         Double_t cti   = ct();
    374                         //
    375                         Double_t Ri    = rh[ii];
    376                         Double_t ArgRp = (Ri*Ci + (1 + Ci * Di)*Di / Ri) / (1 + 2 * Ci*Di);
    377                         Double_t ArgRz = Ci * TMath::Sqrt((Ri*Ri - Di * Di) / (1 + 2 * Ci*Di));
    378                         TVectorD dRphi(5); dRphi.Zero();                // R-phi derivatives @ const. R
    379                         TVectorD dRz(5); dRz.Zero();                    // z     derivatives @ const. R
    380                         //
    381                         // Derivative overflow protection
    382                         Double_t dMin = 0.8;
    383                         dRphi(0) = (1 - 2 * Ci*Ci*Ri*Ri) /
    384                                 TMath::Max(dMin,TMath::Sqrt(1 - ArgRp * ArgRp));        // D derivative
    385                         dRphi(1) = Ri;                                                                                                          // phi0 derivative
    386                         dRphi(2) = Ri * Ri /
    387                                 TMath::Max(dMin,TMath::Sqrt(1 - ArgRp * ArgRp));                                // C derivative
    388                         dRphi(3) = 0.0;                                                                                         // z0 derivative
    389                         dRphi(4) = 0.0;                                                                                         // cot(theta) derivative
    390                         //
    391                         dRz(0) = -cti * Di /
    392                                 (Ri*TMath::Max(dMin,TMath::Sqrt(1 - Ci * Ci*Ri*Ri)));   // D
    393                         dRz(1) = 0.0;                                                                                           // Phi0
    394                         dRz(2) = cti * (Ri*Ci / TMath::Sqrt(1-Ri*Ri*Ci*Ci) -            // C
    395                                 TMath::ASin(Ri*Ci)) / (Ci*Ci);
    396                         dRz(3) = 1.0;                                                                                           // Z0
    397                         dRz(4) = TMath::ASin(ArgRz) / Ci;                                                       // Cot(theta)
    398                         //
    399                         for (Int_t nmi = 0; nmi < nmeai; nmi++)
    400                         {
    401                                 mTl[im] = i;
    402                                 Double_t stri = 0;
    403                                 Double_t sig = 0;
    404                                 if (nmi + 1 == 1)               // Upper layer measurements
    405                                 {
    406                                         stri = fG->lStU(i);     // Stereo angle
    407                                         Double_t csa = TMath::Cos(stri);
    408                                         Double_t ssa = TMath::Sin(stri);
    409                                         sig = fG->lSgU(i);      // Resolution
    410                                         if (ityp == 1)          // Barrel type layer (Measure R-phi, stereo or z at const. R)
    411                                         {
    412                                                 //
    413                                                 // Almost exact solution (CD<<1, D<<R)
    414                                                 Rm(im, 0) = csa * dRphi(0) - ssa * dRz(0);      // D derivative
    415                                                 Rm(im, 1) = csa * dRphi(1) - ssa * dRz(1);      // phi0 derivative
    416                                                 Rm(im, 2) = csa * dRphi(2) - ssa * dRz(2);      // C derivative
    417                                                 Rm(im, 3) = csa * dRphi(3) - ssa * dRz(3);      // z0 derivative
    418                                                 Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4);      // cot(theta) derivative
    419                                         }
    420                                         if (ityp == 2)          // Z type layer (Measure phi at const. Z)
    421                                         {
    422                                                 Rm(im, 0) = 1.0;                                        // D derivative
    423                                                 Rm(im, 1) = rh[ii];                                     // phi0 derivative
    424                                                 Rm(im, 2) = rh[ii] * rh[ii];            // C derivative
    425                                                 Rm(im, 3) = 0;                                          // z0 derivative
    426                                                 Rm(im, 4) = 0;                                          // cot(theta) derivative
    427                                         }
    428                                 }
    429                                 if (nmi + 1 == 2)                       // Lower layer measurements
    430                                 {
    431                                         stri = fG->lStL(i);             // Stereo angle
    432                                         Double_t csa = TMath::Cos(stri);
    433                                         Double_t ssa = TMath::Sin(stri);
    434                                         sig = fG->lSgL(i);              // Resolution
    435                                         if (ityp == 1)                  // Barrel type layer (measure R-phi, stereo or z at const. R)
    436                                         {
    437                                                 //
    438                                                 // Almost exact solution (CD<<1, D<<R)
    439                                                 Rm(im, 0) = csa * dRphi(0) - ssa * dRz(0);      // D derivative
    440                                                 Rm(im, 1) = csa * dRphi(1) - ssa * dRz(1);      // phi0 derivative
    441                                                 Rm(im, 2) = csa * dRphi(2) - ssa * dRz(2);      // C derivative
    442                                                 Rm(im, 3) = csa * dRphi(3) - ssa * dRz(3);      // z0 derivative
    443                                                 Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4);      // cot(theta) derivative
    444                                         }
    445                                         if (ityp == 2)                  // Z type layer (Measure R at const. z)
    446                                         {
    447                                                 Rm(im, 0) = 0;                                                          // D derivative
    448                                                 Rm(im, 1) = 0;                                                          // phi0 derivative
    449                                                 Rm(im, 2) = 0;                                                          // C derivative
    450                                                 Rm(im, 3) = -1.0 / ct();                                        // z0 derivative
    451                                                 Rm(im, 4) = -rh[ii] / ct();                                     // cot(theta) derivative
    452                                         }
    453                                 }
    454                                 // Derivative calculation completed
    455                                 //
    456                                 // Now calculate measurement error matrix
    457                                 //
    458                                 Int_t km = 0;
    459                                 for (Int_t kk = 0; kk <= ii; kk++)
    460                                 {
    461                                         Int_t k = ih[kk];                                       // True layer number
    462                                         Int_t ktyp = fG->lTyp(k);                       // Layer type Barrel or
    463                                         Int_t nmeak = fG->lND(k);                       // # measurements in layer
    464                                         if (fG->isMeasure(k) && nmeak > 0 &&cs[kk] > csMin)
    465                                         {
    466                                                 for (Int_t nmk = 0; nmk < nmeak; nmk++)
    467                                                 {
    468                                                         Double_t strk = 0;
    469                                                         if (nmk + 1 == 1) strk = fG->lStU(k);   // Stereo angle
    470                                                         if (nmk + 1 == 2) strk = fG->lStL(k);   // Stereo angle
    471                                                         if (im == km && Res) Sm(im, km) += sig*sig;     // Detector resolution on diagonal
    472                                                         //
    473                                                         // Loop on all layers below for MS contributions
    474                                                         for (Int_t jj = 0; jj < kk; jj++)
    475                                                         {
    476                                                                 Double_t di = dik(ii, jj);
    477                                                                 Double_t dk = dik(kk, jj);
    478                                                                 Double_t ms = thms[jj];
    479                                                                 Double_t msk = ms; Double_t msi = ms;
    480                                                                 if (ityp == 1) msi = ms / snt;                  // Barrel
    481                                                                 else if (ityp == 2) msi = ms / cst;             // Disk
    482                                                                 if (ktyp == 1) msk = ms / snt;                  // Barrel
    483                                                                 else if (ktyp == 2) msk = ms / cst;             // Disk
    484                                                                 Double_t ci = TMath::Cos(stri); Double_t si = TMath::Sin(stri);
    485                                                                 Double_t ck = TMath::Cos(strk); Double_t sk = TMath::Sin(strk);
    486                                                                 Sm(im, km) += di*dk*(ci*ck*ms*ms + si*sk*msi*msk);                      // Ms contribution
    487                                                         }
    488                                                         //
    489                                                         Sm(km, im) = Sm(im, km);
    490                                                         km++;
    491                                                 }
    492                                         }
    493                                 }
    494                                 im++; mTot = im;
    495                         }
    496                 }
    497         }
    498         Sm.ResizeTo(mTot, mTot);
    499         Rm.ResizeTo(mTot, 5);
    500         //
    501         //**********************************************************************
    502         // Calculate covariance from derivatives and measurement error matrix  *
    503         //**********************************************************************
    504         //
    505         TMatrixDSym DSmInv(mTot); DSmInv.Zero();
    506         for (Int_t id = 0; id < mTot; id++) DSmInv(id, id) = 1.0 / TMath::Sqrt(Sm(id, id));
    507         TMatrixDSym SmN = Sm.Similarity(DSmInv);        // Normalize diagonal to 1
    508         //SmN.Invert();
    509         //
    510         // Protected matrix inversions
    511         //
    512         TDecompChol Chl(SmN);
    513         TMatrixDSym SmNinv = SmN;
    514         if (Chl.Decompose())
    515         {
    516                 Bool_t OK;
    517                 SmNinv = Chl.Invert(OK);
    518         }
    519         else
    520         {
    521                 cout << "SolTrack::CovCalc: Error matrix not positive definite. Recovering ...." << endl;
    522                 if (ntry < ntrymax)
    523                 {
    524                         SmNinv.Print();
    525                         ntry++;
    526                 }
    527                 TMatrixDSym rSmN = MakePosDef(SmN); SmN = rSmN;
    528                 TDecompChol rChl(SmN);
    529                 SmNinv = SmN;
    530                 Bool_t OK = rChl.Decompose();
    531                 SmNinv    = rChl.Invert(OK);
    532         }
    533         Sm = SmNinv.Similarity(DSmInv);                 // Error matrix inverted
    534         TMatrixDSym H = Sm.SimilarityT(Rm);             // Calculate half Hessian
    535         // Normalize before inversion
    536         const Int_t Npar = 5;
    537         TMatrixDSym DHinv(Npar); DHinv.Zero();
    538         for (Int_t i = 0; i < Npar; i++)DHinv(i, i) = 1.0 / TMath::Sqrt(H(i, i));
    539         TMatrixDSym Hnrm = H.Similarity(DHinv);
    540         // Invert and restore
    541         Hnrm.Invert();
    542         fCov = Hnrm.Similarity(DHinv);
     166  // Input flags:
     167  //        Res = .TRUE. turn on resolution effects/Use standard resolutions
     168  //            .FALSE. set all resolutions to 0
     169  //        MS  = .TRUE. include Multiple Scattering
     170  // Assumptions:
     171  // 1. Measurement layers can do one or two measurements
     172  // 2. On disks at constant z:
     173  //    - Upper side measurement is phi
     174  //    - Lower side measurement is R
     175
     176  // Fill list of layers hit
     177  Int_t ntry = 0;
     178  Int_t ntrymax = 0;
     179  Int_t Nhit = nHit(); // Total number of layers hit
     180  Double_t *zhh = new Double_t[Nhit]; // z of hit
     181  Double_t *rhh = new Double_t[Nhit]; // r of hit
     182  Double_t *dhh = new Double_t[Nhit]; // distance of hit from origin
     183  Int_t    *ihh = new Int_t[Nhit];    // true index of layer
     184  Int_t kmh; // Number of measurement layers hit
     185
     186  kmh = HitList(ihh, rhh, zhh); // hit layer list
     187  Int_t mTot = 0; // Total number of measurements
     188  for (Int_t i = 0; i < Nhit; i++)
     189  {
     190    dhh[i] = TMath::Sqrt(rhh[i] * rhh[i] + zhh[i] * zhh[i]);
     191    if (fG->isMeasure(ihh[i])) mTot += fG->lND(ihh[i]); // Count number of measurements
     192  }
     193  // Order hit list by increasing distance from origin
     194  Int_t *hord = new Int_t[Nhit]; // hit order by increasing distance from origin
     195  TMath::Sort(Nhit, dhh, hord, kFALSE); // Order by increasing distance from origin
     196  Double_t *zh = new Double_t[Nhit]; // d-ordered z of hit
     197  Double_t *rh = new Double_t[Nhit]; // d-ordered r of hit
     198  Int_t    *ih = new Int_t[Nhit];    // d-ordered true index of layer
     199  for (Int_t i = 0; i < Nhit; i++)
     200  {
     201    Int_t il = hord[i]; // Hit layer numbering
     202    zh[i] = zhh[il];
     203    rh[i] = rhh[il];
     204    ih[i] = ihh[il];
     205  }
     206  // Store interdistances and multiple scattering angles
     207  Double_t sn2t = 1.0 / (1 + ct()*ct()); //sin^2 theta of track
     208  Double_t cs2t = 1.0 - sn2t;            //cos^2 theta
     209  Double_t snt = TMath::Sqrt(sn2t);      // sin theta
     210  Double_t cst = TMath::Sqrt(cs2t);      // cos theta
     211  Double_t px0 = pt() * TMath::Cos(phi0()); // Momentum at minimum approach
     212  Double_t py0 = pt() * TMath::Sin(phi0());
     213  Double_t pz0 = pt() * ct();
     214  TMatrixDSym dik(Nhit);  dik.Zero(); // Distances between layers
     215  Double_t *thms = new Double_t[Nhit]; // Scattering angles/plane
     216  Double_t *cs = new Double_t[Nhit]; // Cosine of angle with layer normal
     217  for (Int_t ii = 0; ii < Nhit; ii++) // Hit layer loop
     218  {
     219    Int_t i = ih[ii]; // Get true layer number
     220    Double_t B = C()*TMath::Sqrt((rh[ii] * rh[ii] - D()*D()) / (1 + 2 * C()*D()));
     221    Double_t pxi = px0*(1-2*B*B)-2*py0*B*TMath::Sqrt(1-B*B); // Momentum at scattering layer
     222    Double_t pyi = py0*(1-2*B*B)+2*px0*B*TMath::Sqrt(1-B*B);
     223    Double_t pzi = pz0;
     224    Double_t ArgRp = (rh[ii]*C() + (1 + C() * D())*D() / rh[ii]) / (1 + 2 * C()*D());
     225    Double_t phi = phi0() + TMath::ASin(ArgRp);
     226    Double_t nx = TMath::Cos(phi); // Barrel layer normal
     227    Double_t ny = TMath::Sin(phi);
     228    Double_t nz = 0.0;
     229    if (fG->lTyp(i) == 2) // this is Z layer
     230    {
     231      nx = 0.0;
     232      ny = 0.0;
     233      nz = 1.0;
     234    }
     235    Double_t corr = (pxi*nx + pyi * ny + pzi * nz) / p();
     236    cs[ii] = corr;
     237    Double_t Rlf = fG->lTh(i) / (corr*fG->lX0(i)); // Rad. length fraction
     238    thms[ii] = 0.0136*TMath::Sqrt(Rlf)*(1.0 + 0.038*TMath::Log(Rlf)) / p(); // MS angle
     239    if (!MS)thms[ii] = 0;
     240    //
     241    for (Int_t kk = 0; kk < ii; kk++) // Fill distances between layers
     242    {
     243      Double_t Ci = C();
     244      dik(ii, kk) = (TMath::ASin(Ci*rh[ii])-TMath::ASin(Ci*rh[kk]))/(Ci*snt);
     245      dik(kk, ii) = dik(ii, kk);
     246    }
     247    // Store momentum components for resolution correction cosines
     248    Double_t *pRi = new Double_t[Nhit];
     249    pRi[ii] = TMath::Abs(pxi * TMath::Cos(phi) + pyi * TMath::Sin(phi)); // Radial component
     250    Double_t *pPhi = new Double_t[Nhit];
     251    pPhi[ii] = TMath::Abs(pxi * TMath::Sin(phi) - pyi * TMath::Cos(phi)); // Phi component
     252  }
     253  // Fill measurement covariance
     254  Int_t *mTl = new Int_t[mTot]; // Pointer from measurement number to true layer number
     255  TMatrixDSym Sm(mTot); Sm.Zero(); // Measurement covariance
     256  TMatrixD Rm(mTot, 5); // Derivative matrix
     257  Int_t im = 0;
     258  // Fill derivatives and error matrix with MS
     259  Double_t AngMax = 90.; Double_t AngMx = AngMax * TMath::Pi() / 180.;
     260  Double_t csMin = TMath::Cos(AngMx); // Set maximum angle wrt normal
     261  //
     262  for (Int_t ii = 0; ii < Nhit; ii++)
     263  {
     264    Int_t i = ih[ii]; // True layer number
     265    Int_t ityp  = fG->lTyp(i); // Layer type Barrel or Z
     266    Int_t nmeai = fG->lND(i); // # measurements in layer
     267    if (fG->isMeasure(i) && nmeai >0 && cs[ii] > csMin)
     268    {
     269      Double_t Di    = D(); // Get true track parameters
     270      Double_t phi0i = phi0();
     271      Double_t Ci    = C();
     272      Double_t z0i   = z0();
     273      Double_t cti   = ct();
     274      //
     275      Double_t Ri    = rh[ii];
     276      Double_t ArgRp = (Ri*Ci + (1 + Ci * Di)*Di / Ri) / (1 + 2 * Ci*Di);
     277      Double_t ArgRz = Ci * TMath::Sqrt((Ri*Ri - Di * Di) / (1 + 2 * Ci*Di));
     278      TVectorD dRphi(5); dRphi.Zero(); // R-phi derivatives @ const. R
     279      TVectorD dRz(5); dRz.Zero(); // z derivatives @ const. R
     280      // Derivative overflow protection
     281      Double_t dMin = 0.8;
     282      dRphi(0) = (1 - 2 * Ci*Ci*Ri*Ri) /
     283        TMath::Max(dMin,TMath::Sqrt(1 - ArgRp * ArgRp)); // D derivative
     284      dRphi(1) = Ri; // phi0 derivative
     285      dRphi(2) = Ri * Ri /
     286        TMath::Max(dMin,TMath::Sqrt(1 - ArgRp * ArgRp)); // C derivative
     287      dRphi(3) = 0.0; // z0 derivative
     288      dRphi(4) = 0.0; // cot(theta) derivative
     289
     290      dRz(0) = -cti * Di /
     291        (Ri*TMath::Max(dMin,TMath::Sqrt(1 - Ci * Ci*Ri*Ri))); // D
     292      dRz(1) = 0.0; // Phi0
     293      dRz(2) = cti * (Ri*Ci / TMath::Sqrt(1-Ri*Ri*Ci*Ci) -
     294        TMath::ASin(Ri*Ci)) / (Ci*Ci); // C
     295      dRz(3) = 1.0; // Z0
     296      dRz(4) = TMath::ASin(ArgRz) / Ci; // Cot(theta)
     297
     298      for (Int_t nmi = 0; nmi < nmeai; nmi++)
     299      {
     300        mTl[im] = i;
     301        Double_t stri = 0;
     302        Double_t sig = 0;
     303        if (nmi + 1 == 1) // Upper layer measurements
     304        {
     305          stri = fG->lStU(i); // Stereo angle
     306          Double_t csa = TMath::Cos(stri);
     307          Double_t ssa = TMath::Sin(stri);
     308          sig = fG->lSgU(i); // Resolution
     309          if (ityp == 1) // Barrel type layer (Measure R-phi, stereo or z at const. R)
     310          {
     311            // Almost exact solution (CD<<1, D<<R)
     312            Rm(im, 0) = csa * dRphi(0) - ssa * dRz(0); // D derivative
     313            Rm(im, 1) = csa * dRphi(1) - ssa * dRz(1); // phi0 derivative
     314            Rm(im, 2) = csa * dRphi(2) - ssa * dRz(2); // C derivative
     315            Rm(im, 3) = csa * dRphi(3) - ssa * dRz(3); // z0 derivative
     316            Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4); // cot(theta) derivative
     317          }
     318          if (ityp == 2) // Z type layer (Measure phi at const. Z)
     319          {
     320            Rm(im, 0) = 1.0; // D derivative
     321            Rm(im, 1) = rh[ii]; // phi0 derivative
     322            Rm(im, 2) = rh[ii] * rh[ii]; // C derivative
     323            Rm(im, 3) = 0; // z0 derivative
     324            Rm(im, 4) = 0; // cot(theta) derivative
     325          }
     326        }
     327        if (nmi + 1 == 2) // Lower layer measurements
     328        {
     329          stri = fG->lStL(i); // Stereo angle
     330          Double_t csa = TMath::Cos(stri);
     331          Double_t ssa = TMath::Sin(stri);
     332          sig = fG->lSgL(i); // Resolution
     333          if (ityp == 1) // Barrel type layer (measure R-phi, stereo or z at const. R)
     334          {
     335            // Almost exact solution (CD<<1, D<<R)
     336            Rm(im, 0) = csa * dRphi(0) - ssa * dRz(0); // D derivative
     337            Rm(im, 1) = csa * dRphi(1) - ssa * dRz(1); // phi0 derivative
     338            Rm(im, 2) = csa * dRphi(2) - ssa * dRz(2); // C derivative
     339            Rm(im, 3) = csa * dRphi(3) - ssa * dRz(3); // z0 derivative
     340            Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4); // cot(theta) derivative
     341          }
     342          if (ityp == 2) // Z type layer (Measure R at const. z)
     343          {
     344            Rm(im, 0) = 0; // D derivative
     345            Rm(im, 1) = 0; // phi0 derivative
     346            Rm(im, 2) = 0; // C derivative
     347            Rm(im, 3) = -1.0 / ct(); // z0 derivative
     348            Rm(im, 4) = -rh[ii] / ct(); // cot(theta) derivative
     349          }
     350        }
     351        // Derivative calculation completed
     352        // Now calculate measurement error matrix
     353        Int_t km = 0;
     354        for (Int_t kk = 0; kk <= ii; kk++)
     355        {
     356          Int_t k = ih[kk]; // True layer number
     357          Int_t ktyp = fG->lTyp(k); // Layer type Barrel or
     358          Int_t nmeak = fG->lND(k); // # measurements in layer
     359          if (fG->isMeasure(k) && nmeak > 0 &&cs[kk] > csMin)
     360          {
     361            for (Int_t nmk = 0; nmk < nmeak; nmk++)
     362            {
     363              Double_t strk = 0;
     364              if (nmk + 1 == 1) strk = fG->lStU(k); // Stereo angle
     365              if (nmk + 1 == 2) strk = fG->lStL(k); // Stereo angle
     366              if (im == km && Res) Sm(im, km) += sig*sig; // Detector resolution on diagonal
     367              //
     368              // Loop on all layers below for MS contributions
     369              for (Int_t jj = 0; jj < kk; jj++)
     370              {
     371                Double_t di = dik(ii, jj);
     372                Double_t dk = dik(kk, jj);
     373                Double_t ms = thms[jj];
     374                Double_t msk = ms; Double_t msi = ms;
     375                if (ityp == 1) msi = ms / snt; // Barrel
     376                else if (ityp == 2) msi = ms / cst; // Disk
     377                if (ktyp == 1) msk = ms / snt; // Barrel
     378                else if (ktyp == 2) msk = ms / cst; // Disk
     379                Double_t ci = TMath::Cos(stri); Double_t si = TMath::Sin(stri);
     380                Double_t ck = TMath::Cos(strk); Double_t sk = TMath::Sin(strk);
     381                Sm(im, km) += di*dk*(ci*ck*ms*ms + si*sk*msi*msk); // Ms contribution
     382              }
     383              Sm(km, im) = Sm(im, km);
     384              km++;
     385            }
     386          }
     387        }
     388        im++; mTot = im;
     389      }
     390    }
     391  }
     392  Sm.ResizeTo(mTot, mTot);
     393  Rm.ResizeTo(mTot, 5);
     394
     395  // Calculate covariance from derivatives and measurement error matrix
     396  TMatrixDSym DSmInv(mTot); DSmInv.Zero();
     397  for (Int_t id = 0; id < mTot; id++) DSmInv(id, id) = 1.0 / TMath::Sqrt(Sm(id, id));
     398  TMatrixDSym SmN = Sm.Similarity(DSmInv);  // Normalize diagonal to 1
     399  // Protected matrix inversions
     400  TDecompChol Chl(SmN);
     401  TMatrixDSym SmNinv = SmN;
     402  if (Chl.Decompose())
     403  {
     404    Bool_t OK;
     405    SmNinv = Chl.Invert(OK);
     406  }
     407  else
     408  {
     409    cout << "SolTrack::CovCalc: Error matrix not positive definite. Recovering ...." << endl;
     410    if (ntry < ntrymax)
     411    {
     412      SmNinv.Print();
     413      ntry++;
     414    }
     415    TMatrixDSym rSmN = MakePosDef(SmN); SmN = rSmN;
     416    TDecompChol rChl(SmN);
     417    SmNinv = SmN;
     418    Bool_t OK = rChl.Decompose();
     419    SmNinv    = rChl.Invert(OK);
     420  }
     421  Sm = SmNinv.Similarity(DSmInv); // Error matrix inverted
     422  TMatrixDSym H = Sm.SimilarityT(Rm); // Calculate half Hessian
     423  // Normalize before inversion
     424  const Int_t Npar = 5;
     425  TMatrixDSym DHinv(Npar); DHinv.Zero();
     426  for (Int_t i = 0; i < Npar; i++)DHinv(i, i) = 1.0 / TMath::Sqrt(H(i, i));
     427  TMatrixDSym Hnrm = H.Similarity(DHinv);
     428  // Invert and restore
     429  Hnrm.Invert();
     430  fCov = Hnrm.Similarity(DHinv);
    543431}
    544432
    … …  
    546434TMatrixDSym SolTrack::MakePosDef(TMatrixDSym NormMat)
    547435{
    548         //
    549         // Input: symmetric matrix with 1's on diagonal
    550         // Output: positive definite matrix with 1's on diagonal
    551         //
    552         // Default return value
    553         TMatrixDSym rMatN = NormMat;
    554         // Check the diagonal
    555         Bool_t Check = kFALSE;
    556         Int_t Size = NormMat.GetNcols();
    557         for (Int_t i = 0; i < Size; i++)if (TMath::Abs(NormMat(i, i) - 1.0)>1.0E-15)Check = kTRUE;
    558         if (Check)
    559         {
    560                 cout << "SolTrack::MakePosDef: input matrix doesn ot have 1 on diagonal. Abort." << endl;
    561                 return rMatN;
    562         }
    563         //
    564         // Diagonalize matrix
    565         TMatrixDSymEigen Eign(NormMat);
    566         TMatrixD U = Eign.GetEigenVectors();
    567         TVectorD lambda = Eign.GetEigenValues();
    568         // Reset negative eigenvalues to small positive value
    569         TMatrixDSym D(Size); D.Zero(); Double_t eps = 1.0e-13;
    570         for (Int_t i = 0; i < Size; i++)
    571         {
    572                 D(i, i) = lambda(i);
    573                 if (lambda(i) <= 0) D(i, i) = eps;
    574         }
    575          //Rebuild matrix
    576         TMatrixD Ut(TMatrixD::kTransposed, U);
    577         TMatrixD rMat = (U*D)*Ut;                               // Now it is positive defite
    578         // Restore all ones on diagonal
    579         for (Int_t i1 = 0; i1 < Size; i1++)
    580         {
    581                 Double_t rn1 = TMath::Sqrt(rMat(i1, i1));
    582                 for (Int_t i2 = 0; i2 <= i1; i2++)
    583                 {
    584                         Double_t rn2 = TMath::Sqrt(rMat(i2, i2));
    585                         rMatN(i1, i2) = 0.5*(rMat(i1, i2) + rMat(i2, i1)) / (rn1*rn2);
    586                         rMatN(i2, i1) = rMatN(i1, i2);
    587                 }
    588         }
    589         return rMatN;
    590 }
     436  // Input: symmetric matrix with 1's on diagonal
     437  // Output: positive definite matrix with 1's on diagonal
     438
     439  // Default return value
     440  TMatrixDSym rMatN = NormMat;
     441  // Check the diagonal
     442  Bool_t Check = kFALSE;
     443  Int_t Size = NormMat.GetNcols();
     444  for (Int_t i = 0; i < Size; i++)if (TMath::Abs(NormMat(i, i) - 1.0)>1.0E-15)Check = kTRUE;
     445  if (Check)
     446  {
     447    cout << "SolTrack::MakePosDef: input matrix doesn't have 1 on diagonal. Abort." << endl;
     448    return rMatN;
     449  }
     450  // Diagonalize matrix
     451  TMatrixDSymEigen Eign(NormMat);
     452  TMatrixD U = Eign.GetEigenVectors();
     453  TVectorD lambda = Eign.GetEigenValues();
     454  // Reset negative eigenvalues to small positive value
     455  TMatrixDSym D(Size); D.Zero(); Double_t eps = 1.0e-13;
     456  for (Int_t i = 0; i < Size; i++)
     457  {
     458    D(i, i) = lambda(i);
     459    if (lambda(i) <= 0) D(i, i) = eps;
     460  }
     461  // Rebuild matrix
     462  TMatrixD Ut(TMatrixD::kTransposed, U);
     463  TMatrixD rMat = (U*D)*Ut;       // Now it is positive defite
     464  // Restore all ones on diagonal
     465  for (Int_t i1 = 0; i1 < Size; i1++)
     466  {
     467    Double_t rn1 = TMath::Sqrt(rMat(i1, i1));
     468    for (Int_t i2 = 0; i2 <= i1; i2++)
     469    {
     470      Double_t rn2 = TMath::Sqrt(rMat(i2, i2));
     471      rMatN(i1, i2) = 0.5*(rMat(i1, i2) + rMat(i2, i1)) / (rn1*rn2);
     472      rMatN(i2, i1) = rMatN(i1, i2);
     473    }
     474  }
     475  return rMatN;
     476}

  • external/TrackCovariance/SolTrack.h

    r3e4e196 reee976c2  
    33
    44#include <TMath.h>
    5 #include <TVector3.h>
    65#include <TMatrixDSym.h>
    76
    … …  
    2524  // Constructors
    2625  SolTrack(Double_t *x, Double_t *p, SolGeom *G);
    27         SolTrack(TVector3 x, TVector3 p, SolGeom* G);
    2826  SolTrack(Double_t D, Double_t phi0, Double_t C, Double_t z0, Double_t ct, SolGeom *G);
    2927  // Destructor
    … …  
    5957  // Track hit management
    6058  Int_t nHit();
    61   Int_t nmHit();
    6259  Bool_t HitLayer(Int_t Layer, Double_t &R, Double_t &phi, Double_t &zz);
    6360  Int_t HitList(Int_t *&ihh, Double_t *&rhh, Double_t *&zhh);
    64   Int_t HitListXYZ(Int_t *&ihh, Double_t *&Xh, Double_t *&Yh, Double_t *&Zh);
    65 
    6661  // Make normalized matrix positive definite
    6762  TMatrixDSym MakePosDef(TMatrixDSym NormMat);

  • modules/JetFakeParticle.cc

    r3e4e196 reee976c2  
    146146  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
    147147  {
     148    const TLorentzVector &candidatePosition = candidate->Position;
    148149    const TLorentzVector &candidateMomentum = candidate->Momentum;
    149     eta = candidateMomentum.Eta();
    150     phi = candidateMomentum.Phi();
     150    eta = candidatePosition.Eta();
     151    phi = candidatePosition.Phi();
    151152    pt = candidateMomentum.Pt();
    152153    e = candidateMomentum.E();

  • modules/TrackCovariance.cc

    r3e4e196 reee976c2  
    2121 *  Smears track parameters according to appropriate covariance matrix.
    2222 *
    23  *  \authors F. Bedeschi - INFN Pisa
    24 *            P. Demin - UCLouvain, Louvain-la-Neuve
     23 *  \authors P. Demin - UCLouvain, Louvain-la-Neuve
    2524 *           M. Selvaggi - CERN
    26  *
    2725 *
    2826 */
    … …  
    5250
    5351TrackCovariance::TrackCovariance() :
    54   fGeometry(0), fCovariance(0), fAcx(0), fItInputArray(0)
     52  fGeometry(0), fCovariance(0), fItInputArray(0)
    5553{
    5654  fGeometry = new SolGeom();
    … …  
    7270  fBz = GetDouble("Bz", 0.0);
    7371  fGeometry->Read(GetString("DetectorGeometry", ""));
    74   fNMinHits = GetInt("NMinHits", 6);
    7572
    76   // load geometry
    7773  fCovariance->Calc(fGeometry);
    78   fCovariance->SetMinHits(fNMinHits);
    79   // load geometry
    80   fAcx = fCovariance->AccPnt();
    8174
    8275  // import input array
     76
    8377  fInputArray = ImportArray(GetString("InputArray", "TrackMerger/tracks"));
    8478  fItInputArray = fInputArray->MakeIterator();
    … …  
    10397  Double_t mass, p, pt, q, ct;
    10498  Double_t dd0, ddz, dphi, dct, dp, dpt, dC;
    105 
     99 
    106100
    107101  fItInputArray->Reset();
    … …  
    113107    const TLorentzVector &candidateMomentum = candidate->Momentum;
    114108
    115     if ( !fCovariance->IsAccepted(candidateMomentum.Vect()) ) continue;
    116 
    117109    mass = candidateMomentum.M();
    118110
    … …  
    120112
    121113    mother    = candidate;
    122     candidate = static_cast<Candidate*>(candidate->Clone());
     114    candidate = static_cast<Candidate *>(candidate->Clone());
    123115
    124116    candidate->Momentum.SetVectM(track.GetObsP(), mass);
    125 
     117   
    126118    // converting back to mm
    127119    candidate->InitialPosition.SetXYZT(track.GetObsX().X()*1e03,track.GetObsX().Y()*1e03,track.GetObsX().Z()*1e03,candidatePosition.T()*1e03);
    … …  
    138130    candidate->Yd = track.GetObsX().Y()*1e03;
    139131    candidate->Zd = track.GetObsX().Z()*1e03;
    140 
     132   
    141133    candidate->D0       = track.GetObsPar()[0]*1e03;
    142134    candidate->Phi      = track.GetObsPar()[1];
    … …  
    150142    candidate->Charge   = q;
    151143
    152     dd0       = TMath::Sqrt(track.GetCov()(0, 0))*1e03;
    153     ddz       = TMath::Sqrt(track.GetCov()(3, 3))*1e03;
    154     dphi      = TMath::Sqrt(track.GetCov()(1, 1));
    155     dct       = TMath::Sqrt(track.GetCov()(4, 4));
     144    dd0       = TMath::Sqrt(track.GetCov()(0, 0))*1e03; 
     145    ddz       = TMath::Sqrt(track.GetCov()(3, 3))*1e03; 
     146    dphi      = TMath::Sqrt(track.GetCov()(1, 1)); 
     147    dct       = TMath::Sqrt(track.GetCov()(4, 4)); 
    156148    dpt       = 2 * TMath::Sqrt( track.GetCov()(2, 2))*pt*pt / (0.2998*fBz);
    157149    dp        = TMath::Sqrt((1.+ct*ct)*dpt*dpt + 4*pt*pt*ct*ct*dct*dct/(1.+ct*ct)/(1.+ct*ct));
    … …  
    168160    candidate->TrackResolution = dp / p;
    169161
     162
    170163    candidate->AddCandidate(mother);
    171164

  • modules/TrackCovariance.h

    r3e4e196 reee976c2  
    3636class SolGeom;
    3737class SolGridCov;
    38 class AcceptanceClx;
    3938
    4039class TrackCovariance: public DelphesModule
    … …  
    5049private:
    5150  Double_t fBz;
    52   Int_t fNMinHits;
    5351
    5452  SolGeom *fGeometry;
    5553  SolGridCov *fCovariance;
    56 
    57   AcceptanceClx *fAcx;
    5854
    5955  TIterator *fItInputArray; //!
Note: See TracChangeset for help on using the changeset viewer.