Changes in / [eee976c2:3e4e196] in git
- Files:
-
- 3 added
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
Makefile
reee976c2 r3e4e196 635 635 tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf): \ 636 636 external/Hector/H_VerticalQuadrupole.$(SrcSuf) 637 tmp/external/TrackCovariance/AcceptanceClx.$(ObjSuf): \ 638 external/TrackCovariance/AcceptanceClx.$(SrcSuf) 637 639 tmp/external/TrackCovariance/ObsTrk.$(ObjSuf): \ 638 640 external/TrackCovariance/ObsTrk.$(SrcSuf) … … 1145 1147 tmp/external/Hector/H_VerticalKicker.$(ObjSuf) \ 1146 1148 tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf) \ 1149 tmp/external/TrackCovariance/AcceptanceClx.$(ObjSuf) \ 1147 1150 tmp/external/TrackCovariance/ObsTrk.$(ObjSuf) \ 1148 1151 tmp/external/TrackCovariance/SolGeom.$(ObjSuf) \ -
cards/delphes_card_IDEA.tcl
reee976c2 r3e4e196 1 1 #################################################################### 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 5 5 # 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 # 10 10 ####################################### 11 11 # Order of execution of various modules … … 19 19 MuonTrackingEfficiency 20 20 21 ChargedHadronMomentumSmearing 22 ElectronMomentumSmearing 23 MuonMomentumSmearing 24 25 TrackMerger 21 TrackMergerPre 22 TrackSmearing 23 24 TrackMerger 26 25 Calorimeter 27 26 EFlowMerger … … 34 33 ElectronIsolation 35 34 35 MuonFilter 36 36 MuonEfficiency 37 37 MuonIsolation … … 42 42 GenJetFinder 43 43 GenMissingET 44 44 45 45 FastJetFinder 46 46 … … 98 98 } 99 99 100 # (pt <= 0.1) * (0.00) + 100 # (pt <= 0.1) * (0.00) + 101 101 # (abs(eta) <= 3.0) * (pt > 0.1) * (1.00) + 102 102 # (abs(eta) > 3) * (0.00) … … 119 119 (energy < 0.5 && energy >= 0.3) * (abs(eta) <= 3.0) * (0.65) + 120 120 (energy < 0.3) * (abs(eta) <= 3.0) * (0.06) 121 } 121 } 122 122 } 123 123 … … 140 140 } 141 141 142 143 ########################################144 # Momentum resolution for charged tracks145 ########################################146 147 module MomentumSmearing ChargedHadronMomentumSmearing {148 set InputArray ChargedHadronTrackingEfficiency/chargedHadrons149 set OutputArray chargedHadrons150 151 152 # Resolution given in dpT/pT.153 # IDEAdet154 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 electrons161 ###################################162 163 module MomentumSmearing ElectronMomentumSmearing {164 set InputArray ElectronTrackingEfficiency/electrons165 set OutputArray electrons166 167 # Resolution given in dpT/pT.168 # IDEAdet169 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 muons176 ###############################177 178 module MomentumSmearing MuonMomentumSmearing {179 set InputArray MuonTrackingEfficiency/muons180 set OutputArray muons181 182 # Resolution given in dpT/pT.183 # IDEAdet184 set ResolutionFormula {185 (abs(eta) <= 3.0) * sqrt(0.0001145^2 + 0.0002024^2*pt + (pt*2.093e-005)^2)186 }187 }188 189 142 ############## 190 143 # Track merger 191 144 ############## 192 145 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 193 346 module Merger TrackMerger { 194 347 # add InputArray InputArray 195 add InputArray ChargedHadronMomentumSmearing/chargedHadrons 196 add InputArray ElectronMomentumSmearing/electrons 197 add InputArray MuonMomentumSmearing/muons 348 add InputArray TrackSmearing/tracks 198 349 set OutputArray tracks 199 350 } 200 351 201 352 202 ############# 203 # Calorimeter 204 ############# 353 ############# 354 # Calorimeter 355 ############# 205 356 module DualReadoutCalorimeter Calorimeter { 206 357 set ParticleInputArray ParticlePropagator/stableParticles … … 224 375 set pi [expr {acos(-1)}] 225 376 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: 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: 234 385 set PhiBins {} 235 386 for {set i -120} {$i <= 120} {incr i} { 236 387 add PhiBins [expr {$i * $pi/120}] 237 388 } 238 #deta=0.02 units for |eta| <= 0.88 389 #deta=0.02 units for |eta| <= 0.88 239 390 for {set i -44} {$i < 45} {incr i} { 240 391 set eta [expr {$i * 0.02}] … … 242 393 } 243 394 244 #endcaps: 395 #endcaps: 245 396 set PhiBins {} 246 397 for {set i -120} {$i <= 120} {incr i} { 247 398 add PhiBins [expr {$i* $pi/120}] 248 399 } 249 #deta=0.02 units for 0.88 < |eta| <= 3.0 250 #first, from -3.0 to -0.88 400 #deta=0.02 units for 0.88 < |eta| <= 3.0 401 #first, from -3.0 to -0.88 251 402 for {set i 1} {$i <=106} {incr i} { 252 403 set eta [expr {-3.00 + $i * 0.02}] 253 404 add EtaPhiBins $eta $PhiBins 254 405 } 255 #same for 0.88 to 3.0 406 #same for 0.88 to 3.0 256 407 for {set i 1} {$i <=106} {incr i} { 257 408 set eta [expr {0.88 + $i * 0.02}] … … 259 410 } 260 411 261 # default energy fractions {abs(PDG code)} {Fecal Fhcal} 412 # default energy fractions {abs(PDG code)} {Fecal Fhcal} 262 413 add EnergyFraction {0} {0.0 1.0} 263 # energy fractions for e, gamma and pi0 414 # energy fractions for e, gamma and pi0 264 415 add EnergyFraction {11} {1.0 0.0} 265 416 add EnergyFraction {22} {1.0 0.0} 266 417 add EnergyFraction {111} {1.0 0.0} 267 # energy fractions for muon, neutrinos and neutralinos 418 # energy fractions for muon, neutrinos and neutralinos 268 419 add EnergyFraction {12} {0.0 0.0} 269 420 add EnergyFraction {13} {0.0 0.0} … … 275 426 add EnergyFraction {1000035} {0.0 0.0} 276 427 add EnergyFraction {1000045} {0.0 0.0} 277 # energy fractions for K0short and Lambda 428 # energy fractions for K0short and Lambda 278 429 add EnergyFraction {310} {0.3 0.7} 279 430 add EnergyFraction {3122} {0.3 0.7} 280 431 281 432 282 # set ECalResolutionFormula {resolution formula as a function of eta and energy} 433 # set ECalResolutionFormula {resolution formula as a function of eta and energy} 283 434 set ECalResolutionFormula { 284 435 (abs(eta) <= 0.88 ) * sqrt(energy^2*0.01^2 + energy*0.11^2)+ … … 286 437 } 287 438 288 # set HCalResolutionFormula {resolution formula as a function of eta and energy} 439 # set HCalResolutionFormula {resolution formula as a function of eta and energy} 289 440 set HCalResolutionFormula { 290 441 (abs(eta) <= 0.88 ) * sqrt(energy^2*0.01^2 + energy*0.30^2)+ … … 352 503 } 353 504 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 354 519 ##################### 355 520 # Electron efficiency … … 363 528 364 529 # efficiency formula for electrons 365 set EfficiencyFormula { 530 set EfficiencyFormula { 366 531 (energy < 2.0) * (0.000)+ 367 532 (energy >= 2.0) * (abs(eta) <= 0.88) * (0.99) + … … 388 553 } 389 554 555 390 556 ################# 391 557 # Muon efficiency … … 393 559 394 560 module Efficiency MuonEfficiency { 395 set InputArray Muon MomentumSmearing/muons561 set InputArray MuonFilter/muons 396 562 set OutputArray muons 397 563 … … 399 565 400 566 # efficiency formula for muons 401 set EfficiencyFormula { 567 set EfficiencyFormula { 402 568 (energy < 2.0) * (0.000)+ 403 569 (energy >= 2.0) * (abs(eta) <= 0.88) * (0.99) + … … 547 713 548 714 # add EfficiencyFormula {abs(PDG code)} {efficiency formula as a function of eta and pt} 549 715 550 716 # default efficiency formula (misidentification rate) 551 717 add EfficiencyFormula {0} {0.01} … … 603 769 module TreeWriter TreeWriter { 604 770 # add Branch InputArray BranchName BranchClass 605 771 606 772 add Branch Delphes/allParticles Particle GenParticle 607 773 608 774 add Branch TrackMerger/tracks Track Track 609 775 add Branch Calorimeter/towers Tower Tower 610 776 611 777 add Branch Calorimeter/eflowTracks EFlowTrack Track 612 778 add Branch Calorimeter/eflowPhotons EFlowPhoton Tower … … 616 782 add Branch PhotonEfficiency/photons PhotonEff Photon 617 783 add Branch PhotonIsolation/photons PhotonIso Photon 618 784 619 785 add Branch GenJetFinder/jets GenJet Jet 620 786 add Branch GenMissingET/momentum GenMissingET MissingET 621 787 622 788 add Branch UniqueObjectFinder/jets Jet Jet 623 789 add Branch UniqueObjectFinder/electrons Electron Electron 624 790 add Branch UniqueObjectFinder/photons Photon Photon 625 791 add Branch UniqueObjectFinder/muons Muon Muon 626 627 add Branch JetEnergyScale/jets AntiKtJet Jet 628 792 793 add Branch JetEnergyScale/jets AntiKtJet Jet 794 629 795 add Branch MissingET/momentum MissingET MissingET 630 796 add Branch ScalarHT/energy ScalarHT ScalarHT 631 797 } 632 -
external/TrackCovariance/ObsTrk.cc
reee976c2 r3e4e196 44 44 fCovILC = CovToILC(fCov); 45 45 } 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 46 83 // 47 84 // Destructor … … 98 135 // 99 136 TVector3 Xval; 100 Xval(0) = -D*TMath::Sin(phi0); 101 Xval(1) = D*TMath::Cos(phi0); 137 Xval(0) = -D*TMath::Sin(phi0); 138 Xval(1) = D*TMath::Cos(phi0); 102 139 Xval(2) = z0; 103 140 // … … 135 172 // Check ranges 136 173 Double_t minPt = GC->GetMinPt (); 137 if (pt < minPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << std::endl;174 //if (pt < minPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << std::endl; 138 175 Double_t maxPt = GC->GetMaxPt(); 139 if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl;176 //if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl; 140 177 Double_t minAn = GC->GetMinAng(); 141 if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 142 << " is below grid range of " << minAn << std::endl;178 //if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 179 // << " is below grid range of " << minAn << std::endl; 143 180 Double_t maxAn = GC->GetMaxAng(); 144 if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd145 << " is above grid range of " << maxAn << std::endl;181 //if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 182 // << " is above grid range of " << maxAn << std::endl; 146 183 // 147 184 TMatrixDSym Cov = GC->GetCov(pt, angd); … … 179 216 pACTS(1) = 1000 * Par(3); // z0 from m to mm 180 217 pACTS(2) = Par(1); // Phi0 is unchanged 181 pACTS(3) = TMath::ATan (1.0 / Par(4)) + TMath::PiOver2(); // Theta in [0, pi] range218 pACTS(3) = TMath::ATan2(1.0,Par(4)); // Theta in [0, pi] range 182 219 pACTS(4) = Par(2) / (b*TMath::Sqrt(1 + Par(4)*Par(4))); // q/p in GeV 183 220 pACTS(5) = 0.0; // Time: currently undefined … … 247 284 return cILC; 248 285 } 249 250 251 252 -
external/TrackCovariance/ObsTrk.h
reee976c2 r3e4e196 19 19 // Prefix Gen marks variables before resolution smearing 20 20 // 21 private: 21 private: 22 22 Double_t fB; // Solenoid magnetic field 23 23 SolGridCov *fGC; // Covariance matrix grid … … 25 25 Double_t fObsQ; // Observed track charge 26 26 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 28 28 TVector3 fGenP; // Generated track momentum at track origin 29 29 TVector3 fObsP; // Observed track momentum @ track minimum approach … … 43 43 // 44 44 TVectorD ParToACTS(TVectorD Par); // Parameter conversion 45 TMatrixDSym CovToACTS(TMatrixDSym Cov); // Covariance 45 TMatrixDSym CovToACTS(TMatrixDSym Cov); // Covariance 46 46 // 47 47 // Conversion to ILC parametrization … … 55 55 // x(3) track origin, p(3) track momentum at origin, Q charge, B magnetic field in Tesla 56 56 ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC); // Initialize and generate smeared track 57 // Destructor 57 ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC); // Initialize and generate smeared track 58 // Destructor 58 59 ~ObsTrk(); 59 60 // -
external/TrackCovariance/SolGridCov.cc
reee976c2 r3e4e196 3 3 #include <TMath.h> 4 4 #include <TVectorD.h> 5 #include <TVector3.h> 5 6 #include <TMatrixDSym.h> 6 7 #include <TDecompChol.h> … … 36 37 { 37 38 delete[] fCov; 39 delete fAcc; 38 40 } 39 41 … … 58 60 } 59 61 } 60 } 62 63 // Now make acceptance 64 fAcc = new AcceptanceClx(G); 65 } 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 61 107 // Find bin in grid 62 108 Int_t SolGridCov::GetMinIndex(Double_t xval, Int_t N, TVectorD x) -
external/TrackCovariance/SolGridCov.h
reee976c2 r3e4e196 4 4 #include <TVectorD.h> 5 5 #include <TMatrixDSym.h> 6 #include "AcceptanceClx.h" 6 7 7 8 class SolGeom; … … 17 18 TVectorD fAnga; // Array of angle points in degrees 18 19 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 19 22 // Service routines 20 23 Int_t GetMinIndex(Double_t xval, Int_t N, TVectorD x); // Find bin … … 32 35 Double_t GetMaxAng() { return fAnga(fNang - 1); } 33 36 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 34 46 }; 35 47 -
external/TrackCovariance/SolTrack.cc
reee976c2 r3e4e196 45 45 fCov.ResizeTo(5, 5); 46 46 } 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 } 47 82 48 83 SolTrack::SolTrack(Double_t D, Double_t phi0, Double_t C, Double_t z0, Double_t ct, SolGeom *G) … … 69 104 SolTrack::~SolTrack() 70 105 { 106 delete[] & fp; 107 delete[] & fpar; 71 108 fCov.Clear(); 72 109 } … … 132 169 return kh; 133 170 } 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 134 185 // List of layers hit with intersections 135 186 Int_t SolTrack::HitList(Int_t *&ihh, Double_t *&rhh, Double_t *&zhh) … … 161 212 return kmh; 162 213 } 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 // 163 246 // Covariance matrix estimation 247 // 164 248 void SolTrack::CovCalc(Bool_t Res, Bool_t MS) 165 249 { 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); 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); 431 543 } 432 544 … … 434 546 TMatrixDSym SolTrack::MakePosDef(TMatrixDSym NormMat) 435 547 { 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 } 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 } -
external/TrackCovariance/SolTrack.h
reee976c2 r3e4e196 3 3 4 4 #include <TMath.h> 5 #include <TVector3.h> 5 6 #include <TMatrixDSym.h> 6 7 … … 24 25 // Constructors 25 26 SolTrack(Double_t *x, Double_t *p, SolGeom *G); 27 SolTrack(TVector3 x, TVector3 p, SolGeom* G); 26 28 SolTrack(Double_t D, Double_t phi0, Double_t C, Double_t z0, Double_t ct, SolGeom *G); 27 29 // Destructor … … 57 59 // Track hit management 58 60 Int_t nHit(); 61 Int_t nmHit(); 59 62 Bool_t HitLayer(Int_t Layer, Double_t &R, Double_t &phi, Double_t &zz); 60 63 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 61 66 // Make normalized matrix positive definite 62 67 TMatrixDSym MakePosDef(TMatrixDSym NormMat); -
modules/JetFakeParticle.cc
reee976c2 r3e4e196 146 146 while((candidate = static_cast<Candidate *>(fItInputArray->Next()))) 147 147 { 148 const TLorentzVector &candidatePosition = candidate->Position;149 148 const TLorentzVector &candidateMomentum = candidate->Momentum; 150 eta = candidate Position.Eta();151 phi = candidate Position.Phi();149 eta = candidateMomentum.Eta(); 150 phi = candidateMomentum.Phi(); 152 151 pt = candidateMomentum.Pt(); 153 152 e = candidateMomentum.E(); -
modules/TrackCovariance.cc
reee976c2 r3e4e196 21 21 * Smears track parameters according to appropriate covariance matrix. 22 22 * 23 * \authors P. Demin - UCLouvain, Louvain-la-Neuve 23 * \authors F. Bedeschi - INFN Pisa 24 * P. Demin - UCLouvain, Louvain-la-Neuve 24 25 * M. Selvaggi - CERN 26 * 25 27 * 26 28 */ … … 50 52 51 53 TrackCovariance::TrackCovariance() : 52 fGeometry(0), fCovariance(0), f ItInputArray(0)54 fGeometry(0), fCovariance(0), fAcx(0), fItInputArray(0) 53 55 { 54 56 fGeometry = new SolGeom(); … … 70 72 fBz = GetDouble("Bz", 0.0); 71 73 fGeometry->Read(GetString("DetectorGeometry", "")); 74 fNMinHits = GetInt("NMinHits", 6); 72 75 76 // load geometry 73 77 fCovariance->Calc(fGeometry); 78 fCovariance->SetMinHits(fNMinHits); 79 // load geometry 80 fAcx = fCovariance->AccPnt(); 74 81 75 82 // import input array 76 77 83 fInputArray = ImportArray(GetString("InputArray", "TrackMerger/tracks")); 78 84 fItInputArray = fInputArray->MakeIterator(); … … 97 103 Double_t mass, p, pt, q, ct; 98 104 Double_t dd0, ddz, dphi, dct, dp, dpt, dC; 99 105 100 106 101 107 fItInputArray->Reset(); … … 107 113 const TLorentzVector &candidateMomentum = candidate->Momentum; 108 114 115 if ( !fCovariance->IsAccepted(candidateMomentum.Vect()) ) continue; 116 109 117 mass = candidateMomentum.M(); 110 118 … … 112 120 113 121 mother = candidate; 114 candidate = static_cast<Candidate 122 candidate = static_cast<Candidate*>(candidate->Clone()); 115 123 116 124 candidate->Momentum.SetVectM(track.GetObsP(), mass); 117 125 118 126 // converting back to mm 119 127 candidate->InitialPosition.SetXYZT(track.GetObsX().X()*1e03,track.GetObsX().Y()*1e03,track.GetObsX().Z()*1e03,candidatePosition.T()*1e03); … … 130 138 candidate->Yd = track.GetObsX().Y()*1e03; 131 139 candidate->Zd = track.GetObsX().Z()*1e03; 132 140 133 141 candidate->D0 = track.GetObsPar()[0]*1e03; 134 142 candidate->Phi = track.GetObsPar()[1]; … … 142 150 candidate->Charge = q; 143 151 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)); 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)); 148 156 dpt = 2 * TMath::Sqrt( track.GetCov()(2, 2))*pt*pt / (0.2998*fBz); 149 157 dp = TMath::Sqrt((1.+ct*ct)*dpt*dpt + 4*pt*pt*ct*ct*dct*dct/(1.+ct*ct)/(1.+ct*ct)); … … 160 168 candidate->TrackResolution = dp / p; 161 169 162 163 170 candidate->AddCandidate(mother); 164 171 -
modules/TrackCovariance.h
reee976c2 r3e4e196 36 36 class SolGeom; 37 37 class SolGridCov; 38 class AcceptanceClx; 38 39 39 40 class TrackCovariance: public DelphesModule … … 49 50 private: 50 51 Double_t fBz; 52 Int_t fNMinHits; 51 53 52 54 SolGeom *fGeometry; 53 55 SolGridCov *fCovariance; 56 57 AcceptanceClx *fAcx; 54 58 55 59 TIterator *fItInputArray; //!
Note:
See TracChangeset
for help on using the changeset viewer.