Changes in / [66b1770:90975be] in git

Files:

• Makefile (modified) (2 diffs)
• cards/delphes_card_IDEA.tcl (modified) (20 diffs)
• external/TrackCovariance/AcceptanceClx.cc (added)
• external/TrackCovariance/AcceptanceClx.h (added)
• external/TrackCovariance/ObsTrk.cc (modified) (5 diffs)
• external/TrackCovariance/ObsTrk.h (modified) (4 diffs)
• external/TrackCovariance/SolGridCov.cc (modified) (3 diffs)
• external/TrackCovariance/SolGridCov.h (modified) (3 diffs)
• external/TrackCovariance/SolTrack.cc (modified) (5 diffs)
• external/TrackCovariance/SolTrack.h (modified) (3 diffs)
• modules/TrackCovariance.cc (modified) (9 diffs)
• modules/TrackCovariance.h (modified) (2 diffs)

Makefile

@@ -635 +635 @@
 tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf): \
         external/Hector/H_VerticalQuadrupole.$(SrcSuf)
+tmp/external/TrackCovariance/AcceptanceClx.$(ObjSuf): \
+        external/TrackCovariance/AcceptanceClx.$(SrcSuf)
 tmp/external/TrackCovariance/ObsTrk.$(ObjSuf): \
         external/TrackCovariance/ObsTrk.$(SrcSuf)
@@ -1145 +1147 @@
         tmp/external/Hector/H_VerticalKicker.$(ObjSuf) \
         tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf) \
+        tmp/external/TrackCovariance/AcceptanceClx.$(ObjSuf) \
         tmp/external/TrackCovariance/ObsTrk.$(ObjSuf) \
         tmp/external/TrackCovariance/SolGeom.$(ObjSuf) \

cards/delphes_card_IDEA.tcl

@@ -1 +1 @@
 ####################################################################                                l
-# FCC-ee IDEA detector model                                                                                      
-#                                                                                                    
-# Authors: Elisa Fontanesi, Lorenzo Pezzotti, Massimiliano Antonello                                 
+# FCC-ee IDEA detector model
+#
+# Authors: Elisa Fontanesi, Lorenzo Pezzotti, Massimiliano Antonello
 # email: efontane@bo.infn.it,
-#        lorenzo.pezzotti01@universitadipavia.it,                                                    
-#        m.antonello@uninsubria.it,                                                                  
-#####################################################################                                
-#       
+#        lorenzo.pezzotti01@universitadipavia.it,
+#        m.antonello@uninsubria.it,
+#####################################################################
+#
 #######################################
 # Order of execution of various modules
@@ -19 +19 @@
   MuonTrackingEfficiency
 
-  ChargedHadronMomentumSmearing
-  ElectronMomentumSmearing
-  MuonMomentumSmearing
-
-  TrackMerger 
+  TrackMergerPre
+  TrackSmearing
+
+  TrackMerger
   Calorimeter
   EFlowMerger
@@ -34 +33 @@
   ElectronIsolation
 
+  MuonFilter
   MuonEfficiency
   MuonIsolation
@@ -42 +42 @@
   GenJetFinder
   GenMissingET
-  
+
   FastJetFinder
 
@@ -98 +98 @@
 }
 
-#       (pt <= 0.1)                                     * (0.00) + 
+#       (pt <= 0.1)                                     * (0.00) +
 #       (abs(eta) <= 3.0)               * (pt > 0.1)    * (1.00) +
 #       (abs(eta) > 3)                                  * (0.00)
@@ -119 +119 @@
         (energy < 0.5 && energy >= 0.3) * (abs(eta) <= 3.0)    * (0.65) +
         (energy < 0.3) * (abs(eta) <= 3.0)                     * (0.06)
-    } 
+    }
 }
 
@@ -140 +140 @@
 }
 
-
-########################################
-# Momentum resolution for charged tracks
-########################################
-
-module MomentumSmearing ChargedHadronMomentumSmearing {
-    set InputArray ChargedHadronTrackingEfficiency/chargedHadrons
-    set OutputArray chargedHadrons
-
-
-    # Resolution given in dpT/pT.
-    # IDEAdet
-    set ResolutionFormula {
-        (abs(eta) <= 3.0)                   * sqrt(0.0001145^2 + 0.0002024^2*pt + (pt*2.093e-005)^2)
-    }
-}
-
-###################################
-# Momentum resolution for electrons
-###################################
-
-module MomentumSmearing ElectronMomentumSmearing {
-    set InputArray ElectronTrackingEfficiency/electrons
-    set OutputArray electrons
-
-    # Resolution given in dpT/pT.
-    # IDEAdet
-    set ResolutionFormula {
-        (abs(eta) <= 3.0)                   * sqrt(0.0001145^2 + 0.0002024^2*pt + (pt*2.093e-005)^2) 
-   }  
-}
-
-###############################
-# Momentum resolution for muons
-###############################
-
-module MomentumSmearing MuonMomentumSmearing {
-    set InputArray MuonTrackingEfficiency/muons
-    set OutputArray muons
-
-    # Resolution given in dpT/pT.
-    # IDEAdet
-    set ResolutionFormula {
-        (abs(eta) <= 3.0)                   * sqrt(0.0001145^2 + 0.0002024^2*pt + (pt*2.093e-005)^2) 
-    }
-}
-
 ##############
 # Track merger
 ##############
 
+module Merger TrackMergerPre {
+# add InputArray InputArray
+  add InputArray ChargedHadronTrackingEfficiency/chargedHadrons
+  add InputArray ElectronTrackingEfficiency/electrons
+  add InputArray MuonTrackingEfficiency/muons
+  set OutputArray tracks
+}
+
+
+########################################
+# Smearing for charged tracks
+########################################
+
+module TrackCovariance TrackSmearing {
+    set InputArray TrackMergerPre/tracks
+    set OutputArray tracks
+
+    set Bz 2.0
+
+    ## minimum number of hits to accept a track
+    set NMinHits 6
+
+    ## uses https://raw.githubusercontent.com/selvaggi/FastTrackCovariance/master/GeoIDEA_BASE.txt
+    set DetectorGeometry {
+
+
+      # Layer type 1 = R (barrel) or 2 = z (forward/backward)
+      # Layer label
+      # Minimum dimension z for barrel or R for forward
+      # Maximum dimension z for barrel or R for forward
+      # R/z location of layer
+      # Thickness (meters)
+      # Radiation length (meters)
+      # Number of measurements in layers (1D or 2D)
+      # Stereo angle (rad) - 0(pi/2) = axial(z) layer - Upper side
+      # Stereo angle (rad) - 0(pi/2) = axial(z) layer - Lower side
+      # Resolution Upper side (meters) - 0 = no measurement
+      # Resolution Lower side (meters) - 0 = no measurement
+      # measurement flag = T, scattering only = F
+
+
+      # barrel  name       zmin   zmax   r        w (m)      X0        n_meas  th_up (rad) th_down (rad)    reso_up (m)   reso_down (m)  flag
+
+      1        PIPE       -100    100    0.015    0.0012    0.35276    0        0          0                0             0              0
+      1        VTXLOW     -0.12   0.12   0.017    0.00028   0.0937     2        0          1.5708           3e-006        3e-006         1
+      1        VTXLOW     -0.16   0.16   0.023    0.00028   0.0937     2        0          1.5708           3e-006        3e-006         1
+      1        VTXLOW     -0.16   0.16   0.031    0.00028   0.0937     2        0          1.5708           3e-006        3e-006         1
+      1        VTXHIGH    -1      1      0.32     0.00047   0.0937     2        0          1.5708           7e-006        7e-006         1
+      1        VTXHIGH    -1.05   1.05   0.34     0.00047   0.0937     2        0          1.5708           7e-006        7e-006         1
+
+      # endcap  name       rmin   rmax   z        w (m)      X0        n_meas   th_up (rad)  th_down (rad)   reso_up (m)   reso_down (m) flag
+
+      2        VTXDSK      0.141  0.3   -0.92     0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
+      2        VTXDSK      0.138  0.3   -0.9      0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
+      2        VTXDSK      0.065  0.3   -0.42     0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
+      2        VTXDSK      0.062  0.3   -0.4      0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
+      2        VTXDSK      0.062  0.3    0.4      0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
+      2        VTXDSK      0.065  0.3    0.42     0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
+      2        VTXDSK      0.138  0.3    0.9      0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
+      2        VTXDSK      0.141  0.3    0.92     0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
+
+      1 DCHCANI -2.125 2.125 0.345 0.0002 0.237223 0 0 0 0 0 0
+      1 DCH -2 2 0.36 0.0147748 1400 1 0.0203738 0 0.0001 0 1
+      1 DCH -2 2 0.374775 0.0147748 1400 1 -0.0212097 0 0.0001 0 1
+      1 DCH -2 2 0.38955 0.0147748 1400 1 0.0220456 0 0.0001 0 1
+      1 DCH -2 2 0.404324 0.0147748 1400 1 -0.0228814 0 0.0001 0 1
+      1 DCH -2 2 0.419099 0.0147748 1400 1 0.0237172 0 0.0001 0 1
+      1 DCH -2 2 0.433874 0.0147748 1400 1 -0.024553 0 0.0001 0 1
+      1 DCH -2 2 0.448649 0.0147748 1400 1 0.0253888 0 0.0001 0 1
+      1 DCH -2 2 0.463423 0.0147748 1400 1 -0.0262245 0 0.0001 0 1
+      1 DCH -2 2 0.478198 0.0147748 1400 1 0.0270602 0 0.0001 0 1
+      1 DCH -2 2 0.492973 0.0147748 1400 1 -0.0278958 0 0.0001 0 1
+      1 DCH -2 2 0.507748 0.0147748 1400 1 0.0287314 0 0.0001 0 1
+      1 DCH -2 2 0.522523 0.0147748 1400 1 -0.029567 0 0.0001 0 1
+      1 DCH -2 2 0.537297 0.0147748 1400 1 0.0304025 0 0.0001 0 1
+      1 DCH -2 2 0.552072 0.0147748 1400 1 -0.031238 0 0.0001 0 1
+      1 DCH -2 2 0.566847 0.0147748 1400 1 0.0320734 0 0.0001 0 1
+      1 DCH -2 2 0.581622 0.0147748 1400 1 -0.0329088 0 0.0001 0 1
+      1 DCH -2 2 0.596396 0.0147748 1400 1 0.0337442 0 0.0001 0 1
+      1 DCH -2 2 0.611171 0.0147748 1400 1 -0.0345795 0 0.0001 0 1
+      1 DCH -2 2 0.625946 0.0147748 1400 1 0.0354147 0 0.0001 0 1
+      1 DCH -2 2 0.640721 0.0147748 1400 1 -0.0362499 0 0.0001 0 1
+      1 DCH -2 2 0.655495 0.0147748 1400 1 0.0370851 0 0.0001 0 1
+      1 DCH -2 2 0.67027 0.0147748 1400 1 -0.0379202 0 0.0001 0 1
+      1 DCH -2 2 0.685045 0.0147748 1400 1 0.0387552 0 0.0001 0 1
+      1 DCH -2 2 0.69982 0.0147748 1400 1 -0.0395902 0 0.0001 0 1
+      1 DCH -2 2 0.714595 0.0147748 1400 1 0.0404252 0 0.0001 0 1
+      1 DCH -2 2 0.729369 0.0147748 1400 1 -0.04126 0 0.0001 0 1
+      1 DCH -2 2 0.744144 0.0147748 1400 1 0.0420949 0 0.0001 0 1
+      1 DCH -2 2 0.758919 0.0147748 1400 1 -0.0429296 0 0.0001 0 1
+      1 DCH -2 2 0.773694 0.0147748 1400 1 0.0437643 0 0.0001 0 1
+      1 DCH -2 2 0.788468 0.0147748 1400 1 -0.044599 0 0.0001 0 1
+      1 DCH -2 2 0.803243 0.0147748 1400 1 0.0454336 0 0.0001 0 1
+      1 DCH -2 2 0.818018 0.0147748 1400 1 -0.0462681 0 0.0001 0 1
+      1 DCH -2 2 0.832793 0.0147748 1400 1 0.0471025 0 0.0001 0 1
+      1 DCH -2 2 0.847568 0.0147748 1400 1 -0.0479369 0 0.0001 0 1
+      1 DCH -2 2 0.862342 0.0147748 1400 1 0.0487713 0 0.0001 0 1
+      1 DCH -2 2 0.877117 0.0147748 1400 1 -0.0496055 0 0.0001 0 1
+      1 DCH -2 2 0.891892 0.0147748 1400 1 0.0504397 0 0.0001 0 1
+      1 DCH -2 2 0.906667 0.0147748 1400 1 -0.0512738 0 0.0001 0 1
+      1 DCH -2 2 0.921441 0.0147748 1400 1 0.0521079 0 0.0001 0 1
+      1 DCH -2 2 0.936216 0.0147748 1400 1 -0.0529418 0 0.0001 0 1
+      1 DCH -2 2 0.950991 0.0147748 1400 1 0.0537757 0 0.0001 0 1
+      1 DCH -2 2 0.965766 0.0147748 1400 1 -0.0546095 0 0.0001 0 1
+      1 DCH -2 2 0.980541 0.0147748 1400 1 0.0554433 0 0.0001 0 1
+      1 DCH -2 2 0.995315 0.0147748 1400 1 -0.056277 0 0.0001 0 1
+      1 DCH -2 2 1.01009 0.0147748 1400 1 0.0571106 0 0.0001 0 1
+      1 DCH -2 2 1.02486 0.0147748 1400 1 -0.0579441 0 0.0001 0 1
+      1 DCH -2 2 1.03964 0.0147748 1400 1 0.0587775 0 0.0001 0 1
+      1 DCH -2 2 1.05441 0.0147748 1400 1 -0.0596108 0 0.0001 0 1
+      1 DCH -2 2 1.06919 0.0147748 1400 1 0.0604441 0 0.0001 0 1
+      1 DCH -2 2 1.08396 0.0147748 1400 1 -0.0612773 0 0.0001 0 1
+      1 DCH -2 2 1.09874 0.0147748 1400 1 0.0621104 0 0.0001 0 1
+      1 DCH -2 2 1.11351 0.0147748 1400 1 -0.0629434 0 0.0001 0 1
+      1 DCH -2 2 1.12829 0.0147748 1400 1 0.0637763 0 0.0001 0 1
+      1 DCH -2 2 1.14306 0.0147748 1400 1 -0.0646092 0 0.0001 0 1
+      1 DCH -2 2 1.15784 0.0147748 1400 1 0.0654419 0 0.0001 0 1
+      1 DCH -2 2 1.17261 0.0147748 1400 1 -0.0662746 0 0.0001 0 1
+      1 DCH -2 2 1.18739 0.0147748 1400 1 0.0671071 0 0.0001 0 1
+      1 DCH -2 2 1.20216 0.0147748 1400 1 -0.0679396 0 0.0001 0 1
+      1 DCH -2 2 1.21694 0.0147748 1400 1 0.068772 0 0.0001 0 1
+      1 DCH -2 2 1.23171 0.0147748 1400 1 -0.0696042 0 0.0001 0 1
+      1 DCH -2 2 1.24649 0.0147748 1400 1 0.0704364 0 0.0001 0 1
+      1 DCH -2 2 1.26126 0.0147748 1400 1 -0.0712685 0 0.0001 0 1
+      1 DCH -2 2 1.27604 0.0147748 1400 1 0.0721005 0 0.0001 0 1
+      1 DCH -2 2 1.29081 0.0147748 1400 1 -0.0729324 0 0.0001 0 1
+      1 DCH -2 2 1.30559 0.0147748 1400 1 0.0737642 0 0.0001 0 1
+      1 DCH -2 2 1.32036 0.0147748 1400 1 -0.0745958 0 0.0001 0 1
+      1 DCH -2 2 1.33514 0.0147748 1400 1 0.0754274 0 0.0001 0 1
+      1 DCH -2 2 1.34991 0.0147748 1400 1 -0.0762589 0 0.0001 0 1
+      1 DCH -2 2 1.36468 0.0147748 1400 1 0.0770903 0 0.0001 0 1
+      1 DCH -2 2 1.37946 0.0147748 1400 1 -0.0779215 0 0.0001 0 1
+      1 DCH -2 2 1.39423 0.0147748 1400 1 0.0787527 0 0.0001 0 1
+      1 DCH -2 2 1.40901 0.0147748 1400 1 -0.0795837 0 0.0001 0 1
+      1 DCH -2 2 1.42378 0.0147748 1400 1 0.0804147 0 0.0001 0 1
+      1 DCH -2 2 1.43856 0.0147748 1400 1 -0.0812455 0 0.0001 0 1
+      1 DCH -2 2 1.45333 0.0147748 1400 1 0.0820762 0 0.0001 0 1
+      1 DCH -2 2 1.46811 0.0147748 1400 1 -0.0829068 0 0.0001 0 1
+      1 DCH -2 2 1.48288 0.0147748 1400 1 0.0837373 0 0.0001 0 1
+      1 DCH -2 2 1.49766 0.0147748 1400 1 -0.0845677 0 0.0001 0 1
+      1 DCH -2 2 1.51243 0.0147748 1400 1 0.0853979 0 0.0001 0 1
+      1 DCH -2 2 1.52721 0.0147748 1400 1 -0.086228 0 0.0001 0 1
+      1 DCH -2 2 1.54198 0.0147748 1400 1 0.087058 0 0.0001 0 1
+      1 DCH -2 2 1.55676 0.0147748 1400 1 -0.0878879 0 0.0001 0 1
+      1 DCH -2 2 1.57153 0.0147748 1400 1 0.0887177 0 0.0001 0 1
+      1 DCH -2 2 1.58631 0.0147748 1400 1 -0.0895474 0 0.0001 0 1
+      1 DCH -2 2 1.60108 0.0147748 1400 1 0.0903769 0 0.0001 0 1
+      1 DCH -2 2 1.61586 0.0147748 1400 1 -0.0912063 0 0.0001 0 1
+      1 DCH -2 2 1.63063 0.0147748 1400 1 0.0920356 0 0.0001 0 1
+      1 DCH -2 2 1.64541 0.0147748 1400 1 -0.0928647 0 0.0001 0 1
+      1 DCH -2 2 1.66018 0.0147748 1400 1 0.0936937 0 0.0001 0 1
+      1 DCH -2 2 1.67495 0.0147748 1400 1 -0.0945226 0 0.0001 0 1
+      1 DCH -2 2 1.68973 0.0147748 1400 1 0.0953514 0 0.0001 0 1
+      1 DCH -2 2 1.7045 0.0147748 1400 1 -0.09618 0 0.0001 0 1
+      1 DCH -2 2 1.71928 0.0147748 1400 1 0.0970085 0 0.0001 0 1
+      1 DCH -2 2 1.73405 0.0147748 1400 1 -0.0978369 0 0.0001 0 1
+      1 DCH -2 2 1.74883 0.0147748 1400 1 0.0986651 0 0.0001 0 1
+      1 DCH -2 2 1.7636 0.0147748 1400 1 -0.0994932 0 0.0001 0 1
+      1 DCH -2 2 1.77838 0.0147748 1400 1 0.100321 0 0.0001 0 1
+      1 DCH -2 2 1.79315 0.0147748 1400 1 -0.101149 0 0.0001 0 1
+      1 DCH -2 2 1.80793 0.0147748 1400 1 0.101977 0 0.0001 0 1
+      1 DCH -2 2 1.8227 0.0147748 1400 1 -0.102804 0 0.0001 0 1
+      1 DCH -2 2 1.83748 0.0147748 1400 1 0.103632 0 0.0001 0 1
+      1 DCH -2 2 1.85225 0.0147748 1400 1 -0.104459 0 0.0001 0 1
+      1 DCH -2 2 1.86703 0.0147748 1400 1 0.105286 0 0.0001 0 1
+      1 DCH -2 2 1.8818 0.0147748 1400 1 -0.106113 0 0.0001 0 1
+      1 DCH -2 2 1.89658 0.0147748 1400 1 0.10694 0 0.0001 0 1
+      1 DCH -2 2 1.91135 0.0147748 1400 1 -0.107766 0 0.0001 0 1
+      1 DCH -2 2 1.92613 0.0147748 1400 1 0.108593 0 0.0001 0 1
+      1 DCH -2 2 1.9409 0.0147748 1400 1 -0.109419 0 0.0001 0 1
+      1 DCH -2 2 1.95568 0.0147748 1400 1 0.110246 0 0.0001 0 1
+      1 DCH -2 2 1.97045 0.0147748 1400 1 -0.111072 0 0.0001 0 1
+      1 DCH -2 2 1.98523 0.0147748 1400 1 0.111898 0 0.0001 0 1
+      1 DCH -2 2 2 0.0147748 1400 1 -0.112723 0 0.0001 0 1
+      1 DCHCANO -2.125 2.125 2.02 0.02 1.667 0 0 0 0 0 0
+      1 BSILWRP -2.35 2.35 2.04 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
+      1 BSILWRP -2.35 2.35 2.06 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
+      1 MAG -2.5 2.5 2.25 0.05 0.0658 0 0 0 0 0 0
+      1 BPRESH -2.55 2.55 2.45 0.02 1 2 0 1.5708 7e-005 0.01 1
+
+
+      2 DCHWALL 0.345 2.02 2.125 0.25 5.55 0 0 0 0 0 0
+      2 DCHWALL 0.345 2.02 -2.125 0.25 5.55 0 0 0 0 0 0
+      2 FSILWRP 0.354 2.02 -2.32 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
+      2 FSILWRP 0.35 2.02 -2.3 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
+      2 FSILWRP 0.35 2.02 2.3 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
+      2 FSILWRP 0.354 2.02 2.32 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
+      2 FRAD 0.38 2.09 2.49 0.0043 0.005612 0 0 0 0 0 0
+      2 FRAD 0.38 2.09 -2.49 0.0043 0.005612 0 0 0 0 0 0
+      2 FPRESH 0.39 2.43 -2.55 0.02 1 2 0 1.5708 7e-005 0.01 1
+      2 FPRESH 0.39 2.43 2.55 0.02 1 2 0 1.5708 7e-005 0.01 1
+    }
+
+}
+
+
+##############
+# Track merger
+##############
+
 module Merger TrackMerger {
 # add InputArray InputArray
-  add InputArray ChargedHadronMomentumSmearing/chargedHadrons
-  add InputArray ElectronMomentumSmearing/electrons
-  add InputArray MuonMomentumSmearing/muons
+  add InputArray TrackSmearing/tracks
   set OutputArray tracks
 }
 
 
-#############                                                                                                                         
-# Calorimeter                                                                                                                                            
-#############                                                                                                                                           
+#############
+# Calorimeter
+#############
 module DualReadoutCalorimeter Calorimeter {
   set ParticleInputArray ParticlePropagator/stableParticles
@@ -224 +375 @@
     set pi [expr {acos(-1)}]
 
-    # Lists of the edges of each tower in eta and phi;                                                                                            
-    # each list starts with the lower edge of the first tower;                                                                                                    
-    # the list ends with the higher edged of the last tower.                                                                                        
-    # Barrel:  deta=0.02 towers up to |eta| <= 0.88 ( up to 45°)                                                                                        
-    # Endcaps: deta=0.02 towers up to |eta| <= 3.0 (8.6° = 100 mrad)                                                                                           
-    # Cell size: about 6 cm x 6 cm                                                  
-
-    #barrel:                                                                                        
+    # Lists of the edges of each tower in eta and phi;
+    # each list starts with the lower edge of the first tower;
+    # the list ends with the higher edged of the last tower.
+    # Barrel:  deta=0.02 towers up to |eta| <= 0.88 ( up to 45°)
+    # Endcaps: deta=0.02 towers up to |eta| <= 3.0 (8.6° = 100 mrad)
+    # Cell size: about 6 cm x 6 cm
+
+    #barrel:
     set PhiBins {}
     for {set i -120} {$i <= 120} {incr i} {
         add PhiBins [expr {$i * $pi/120}]
     }
-    #deta=0.02 units for |eta| <= 0.88                                                              
+    #deta=0.02 units for |eta| <= 0.88
     for {set i -44} {$i < 45} {incr i} {
         set eta [expr {$i * 0.02}]
@@ -242 +393 @@
     }
 
-    #endcaps:                                                                                       
+    #endcaps:
     set PhiBins {}
     for {set i -120} {$i <= 120} {incr i} {
         add PhiBins [expr {$i* $pi/120}]
     }
-    #deta=0.02 units for 0.88 < |eta| <= 3.0                                                        
-    #first, from -3.0 to -0.88                                                                      
+    #deta=0.02 units for 0.88 < |eta| <= 3.0
+    #first, from -3.0 to -0.88
     for {set i 1} {$i <=106} {incr i} {
         set eta [expr {-3.00 + $i * 0.02}]
         add EtaPhiBins $eta $PhiBins
     }
-    #same for 0.88 to 3.0                                                                           
+    #same for 0.88 to 3.0
     for  {set i 1} {$i <=106} {incr i} {
         set eta [expr {0.88 + $i * 0.02}]
@@ -259 +410 @@
     }
 
-    # default energy fractions {abs(PDG code)} {Fecal Fhcal}                                                                                  
+    # default energy fractions {abs(PDG code)} {Fecal Fhcal}
     add EnergyFraction {0} {0.0 1.0}
-    # energy fractions for e, gamma and pi0                                                                                                           
+    # energy fractions for e, gamma and pi0
     add EnergyFraction {11} {1.0 0.0}
     add EnergyFraction {22} {1.0 0.0}
     add EnergyFraction {111} {1.0 0.0}
-    # energy fractions for muon, neutrinos and neutralinos                                                                             
+    # energy fractions for muon, neutrinos and neutralinos
     add EnergyFraction {12} {0.0 0.0}
     add EnergyFraction {13} {0.0 0.0}
@@ -275 +426 @@
     add EnergyFraction {1000035} {0.0 0.0}
     add EnergyFraction {1000045} {0.0 0.0}
-    # energy fractions for K0short and Lambda                                                                                                            
+    # energy fractions for K0short and Lambda
     add EnergyFraction {310} {0.3 0.7}
     add EnergyFraction {3122} {0.3 0.7}
 
 
-    # set ECalResolutionFormula {resolution formula as a function of eta and energy}                                
+    # set ECalResolutionFormula {resolution formula as a function of eta and energy}
     set ECalResolutionFormula {
     (abs(eta) <= 0.88 )                     * sqrt(energy^2*0.01^2 + energy*0.11^2)+
@@ -286 +437 @@
     }
 
-    # set HCalResolutionFormula {resolution formula as a function of eta and energy}                                                
+    # set HCalResolutionFormula {resolution formula as a function of eta and energy}
     set HCalResolutionFormula {
     (abs(eta) <= 0.88 )                     * sqrt(energy^2*0.01^2 + energy*0.30^2)+
@@ -352 +503 @@
 }
 
+
+#################
+# Muon filter
+#################
+
+module PdgCodeFilter MuonFilter {
+  set InputArray Calorimeter/eflowTracks
+  set OutputArray muons
+  set Invert true
+  add PdgCode {13}
+  add PdgCode {-13}
+}
+
+
 #####################
 # Electron efficiency
@@ -363 +528 @@
 
   # efficiency formula for electrons
-  set EfficiencyFormula {          
+  set EfficiencyFormula {
         (energy < 2.0)                                         * (0.000)+
         (energy >= 2.0) * (abs(eta) <= 0.88)                   * (0.99) +
@@ -388 +553 @@
 }
 
+
 #################
 # Muon efficiency
@@ -393 +559 @@
 
 module Efficiency MuonEfficiency {
-  set InputArray MuonMomentumSmearing/muons
+  set InputArray MuonFilter/muons
   set OutputArray muons
 
@@ -399 +565 @@
 
   # efficiency formula for muons
-  set EfficiencyFormula {                                   
+  set EfficiencyFormula {
         (energy < 2.0)                                         * (0.000)+
         (energy >= 2.0) * (abs(eta) <= 0.88)                   * (0.99) +
@@ -547 +713 @@
 
   # add EfficiencyFormula {abs(PDG code)} {efficiency formula as a function of eta and pt}
-  
+
   # default efficiency formula (misidentification rate)
   add EfficiencyFormula {0} {0.01}
@@ -603 +769 @@
 module TreeWriter TreeWriter {
     # add Branch InputArray BranchName BranchClass
-    
+
     add Branch Delphes/allParticles Particle GenParticle
 
     add Branch TrackMerger/tracks Track Track
     add Branch Calorimeter/towers Tower Tower
-    
+
     add Branch Calorimeter/eflowTracks EFlowTrack Track
     add Branch Calorimeter/eflowPhotons EFlowPhoton Tower
@@ -616 +782 @@
     add Branch PhotonEfficiency/photons PhotonEff Photon
     add Branch PhotonIsolation/photons PhotonIso Photon
-    
+
     add Branch GenJetFinder/jets GenJet Jet
     add Branch GenMissingET/momentum GenMissingET MissingET
-    
+
     add Branch UniqueObjectFinder/jets Jet Jet
     add Branch UniqueObjectFinder/electrons Electron Electron
     add Branch UniqueObjectFinder/photons Photon Photon
     add Branch UniqueObjectFinder/muons Muon Muon
-    
-    add Branch JetEnergyScale/jets AntiKtJet Jet  
-    
+
+    add Branch JetEnergyScale/jets AntiKtJet Jet
+
     add Branch MissingET/momentum MissingET MissingET
     add Branch ScalarHT/energy ScalarHT ScalarHT
 }
-

external/TrackCovariance/ObsTrk.cc

@@ -44 +44 @@
         fCovILC = CovToILC(fCov);
 }
+
+// x[3] track origin, p[3] track momentum at origin, Q charge, B magnetic field in Tesla
+ObsTrk::ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC)
+{
+        fGC = GC;
+        fGenX.SetXYZ(x[0],x[1],x[2]);
+        fGenP.SetXYZ(p[0],p[1],p[2]);
+        fGenQ = Q;
+        fB = B;
+        fGenPar.ResizeTo(5);
+        fGenParACTS.ResizeTo(6);
+        fGenParILC.ResizeTo(5);
+        fObsPar.ResizeTo(5);
+        fObsParACTS.ResizeTo(6);
+        fObsParILC.ResizeTo(5);
+        fCov.ResizeTo(5, 5);
+        fCovACTS.ResizeTo(6, 6);
+        fCovILC.ResizeTo(5, 5);
+        fGenPar = XPtoPar(fGenX, fGenP, Q);
+        fGenParACTS = ParToACTS(fGenPar);
+        fGenParILC = ParToILC(fGenPar);
+        /*
+        cout << "ObsTrk::ObsTrk: fGenPar";
+        for (Int_t i = 0; i < 5; i++)cout << fGenPar(i) << ", ";
+        cout << endl;
+        */
+        fObsPar = GenToObsPar(fGenPar, fGC);
+        fObsParACTS = ParToACTS(fObsPar);
+        fObsParILC = ParToILC(fObsPar);
+        fObsX = ParToX(fObsPar);
+        fObsP = ParToP(fObsPar);
+        fObsQ = ParToQ(fObsPar);
+        fCovACTS = CovToACTS(fCov);
+        fCovILC = CovToILC(fCov);
+}
+
+
 //
 // Destructor
@@ -98 +135 @@
         //
         TVector3 Xval;
-        Xval(0) = -D*TMath::Sin(phi0); 
-        Xval(1) =  D*TMath::Cos(phi0);  
+        Xval(0) = -D*TMath::Sin(phi0);
+        Xval(1) =  D*TMath::Cos(phi0);
         Xval(2) =  z0;
         //
@@ -135 +172 @@
         // Check ranges
         Double_t minPt = GC->GetMinPt ();
-        if (pt < minPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << std::endl;
+        //if (pt < minPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << std::endl;
         Double_t maxPt = GC->GetMaxPt();
-        if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl;
+        //if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl;
         Double_t minAn = GC->GetMinAng();
-        if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 
-                << " is below grid range of " << minAn << std::endl;
+  //if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
+        //      << " is below grid range of " << minAn << std::endl;
         Double_t maxAn = GC->GetMaxAng();
-        if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
-                << " is above grid range of " << maxAn << std::endl;
+        //if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
+        //      << " is above grid range of " << maxAn << std::endl;
         //
         TMatrixDSym Cov = GC->GetCov(pt, angd);
@@ -179 +216 @@
         pACTS(1) = 1000 * Par(3);       // z0 from m to mm
         pACTS(2) = Par(1);                      // Phi0 is unchanged
-        pACTS(3) = TMath::ATan(1.0 / Par(4)) + TMath::PiOver2();                // Theta in [0, pi] range
+        pACTS(3) = TMath::ATan2(1.0,Par(4));            // Theta in [0, pi] range
         pACTS(4) = Par(2) / (b*TMath::Sqrt(1 + Par(4)*Par(4)));         // q/p in GeV
         pACTS(5) = 0.0;                         // Time: currently undefined
@@ -247 +284 @@
         return cILC;
 }
-
-
-
-        

external/TrackCovariance/ObsTrk.h

@@ -19 +19 @@
         // Prefix Gen marks variables before resolution smearing
         //
-private:        
+private:
         Double_t fB;                                            // Solenoid magnetic field
         SolGridCov *fGC;                                        // Covariance matrix grid
@@ -25 +25 @@
         Double_t fObsQ;                                 // Observed  track charge
         TVector3 fGenX;                                 // Generated track origin (x,y,z)
-        TVector3 fObsX;                                 // Observed  track origin (x,y,z) @ track min. approach 
+        TVector3 fObsX;                                 // Observed  track origin (x,y,z) @ track min. approach
         TVector3 fGenP;                                 // Generated track momentum at track origin
         TVector3 fObsP;                                 // Observed  track momentum @ track minimum approach
@@ -43 +43 @@
         //
         TVectorD ParToACTS(TVectorD Par);               // Parameter conversion
-        TMatrixDSym CovToACTS(TMatrixDSym Cov); // Covariance 
+        TMatrixDSym CovToACTS(TMatrixDSym Cov); // Covariance
         //
         // Conversion to ILC parametrization
@@ -55 +55 @@
         // x(3) track origin, p(3) track momentum at origin, Q charge, B magnetic field in Tesla
         ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC); // Initialize and generate smeared track
-        // Destructor
+        ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC);       // Initialize and generate smeared track
+  // Destructor
         ~ObsTrk();
         //

external/TrackCovariance/SolGridCov.cc

@@ -3 +3 @@
 #include <TMath.h>
 #include <TVectorD.h>
+#include <TVector3.h>
 #include <TMatrixDSym.h>
 #include <TDecompChol.h>
@@ -36 +37 @@
 {
   delete[] fCov;
+  delete fAcc;
 }
 
@@ -58 +60 @@
     }
   }
-}
+
+// Now make acceptance
+fAcc = new AcceptanceClx(G);
+}
+
+
+//
+Bool_t SolGridCov::IsAccepted(Double_t pt, Double_t Theta)
+{
+        //
+        // pt in GeV, Theta in degrees
+        //
+        Bool_t Accept = kFALSE;
+        if (fAcc->HitNumber(pt, Theta) >= fNminHits)Accept = kTRUE;
+        //
+        return Accept;
+}
+//
+Bool_t SolGridCov::IsAccepted(Double_t *p)
+{
+        //
+        // pt in GeV, Theta in degrees
+        //
+        Bool_t Accept = kFALSE;
+        Double_t pt = TMath::Sqrt(p[0] * p[0] + p[1] * p[1]);
+        Double_t th = 180. * TMath::ATan2(pt, p[2])/TMath::Pi();
+        if (fAcc->HitNumber(pt,th) >= fNminHits)Accept = kTRUE;
+        //
+        return Accept;
+}
+//
+Bool_t SolGridCov::IsAccepted(TVector3 p)
+{
+        //
+        // pt in GeV, Theta in degrees
+        //
+        Bool_t Accept = kFALSE;
+        Double_t pt = p.Pt();
+        Double_t th = 180.*TMath::ACos(p.CosTheta())/TMath::Pi();
+        if (fAcc->HitNumber(pt,th) >= fNminHits)Accept = kTRUE;
+        //
+        return Accept;
+}
+
+
 // Find bin in grid
 Int_t SolGridCov::GetMinIndex(Double_t xval, Int_t N, TVectorD x)

external/TrackCovariance/SolGridCov.h

@@ -4 +4 @@
 #include <TVectorD.h>
 #include <TMatrixDSym.h>
+#include "AcceptanceClx.h"
 
 class SolGeom;
@@ -17 +18 @@
   TVectorD fAnga;    // Array of angle points in degrees
   TMatrixDSym *fCov; // Pointers to grid of covariance matrices
+        AcceptanceClx *fAcc;                            // Pointer to acceptance class
+        Int_t fNminHits;                                        // Minimum number of hits to accept track
   // Service routines
   Int_t GetMinIndex(Double_t xval, Int_t N, TVectorD x); // Find bin
@@ -32 +35 @@
   Double_t GetMaxAng() { return fAnga(fNang - 1); }
   TMatrixDSym GetCov(Double_t pt, Double_t ang);
+
+        // Acceptance related methods
+        AcceptanceClx* AccPnt() { return fAcc; };                       // Return Acceptance class pointer
+        void SetMinHits(Int_t MinHits) { fNminHits = MinHits; };        // Set minimum number of hits to accept (default = 6)
+        Int_t GetMinHits() { return fNminHits; };
+        Bool_t IsAccepted(Double_t pt, Double_t Theta);         // From pt (GeV) and theta (degrees)
+        Bool_t IsAccepted(Double_t *p);                                         // From momentum components (GeV)
+        Bool_t IsAccepted(TVector3 p);                                          // As above in Vector3 format
+
 };
 

external/TrackCovariance/SolTrack.cc

@@ -45 +45 @@
   fCov.ResizeTo(5, 5);
 }
+SolTrack::SolTrack(TVector3 x, TVector3 p, SolGeom* G)
+{
+        fG = G;
+        // Store momentum
+        fp[0] = p(0); fp[1] = p(1); fp[2] = p(2);
+        Double_t px = p(0); Double_t py = p(1); Double_t pz = p(2);
+        fx[0] = x(0); fx[1] = x(1); fx[2] = x(2);
+        Double_t xx = x(0); Double_t yy = x(1); Double_t zz = x(2);
+        // Store parameters
+        Double_t pt = TMath::Sqrt(px * px + py * py);
+        Double_t Charge = 1.0;                                          // Don't worry about charge for now
+        Double_t a = -Charge * G->B() * 0.2998;                 // Normalized speed of light
+        Double_t C = a / (2 * pt);                                      // pt in GeV, B in Tesla, C in meters
+        Double_t r2 = xx * xx + yy * yy;
+        Double_t cross = xx * py - yy * px;
+        Double_t T = TMath::Sqrt(pt * pt - 2 * a * cross + a * a * r2);
+        Double_t phi0 = TMath::ATan2((py - a * xx) / T, (px + a * yy) / T);
+        Double_t D;
+        if (pt < 10.0) D = (T - pt) / a;
+        else D = (-2 * cross + a * r2) / (T + pt);
+        Double_t B = C * TMath::Sqrt((r2 - D * D) / (1 + 2 * C * D));
+        Double_t st = TMath::ASin(B) / C;
+        Double_t ct = pz / pt;
+        Double_t z0 = zz - ct * st;
+        fpar[0] = D;
+        fpar[1] = phi0;
+        fpar[2] = C;
+        fpar[3] = z0;
+        fpar[4] = ct;
+        //cout << "SolTrack:: C = " << C << ", fpar[2] = " << fpar[2] << endl;
+        //
+        // Init covariances
+        //
+        fCov.ResizeTo(5, 5);
+}
 
 SolTrack::SolTrack(Double_t D, Double_t phi0, Double_t C, Double_t z0, Double_t ct, SolGeom *G)
@@ -69 +104 @@
 SolTrack::~SolTrack()
 {
+        delete[] & fp;
+        delete[] & fpar;
   fCov.Clear();
 }
@@ -132 +169 @@
   return kh;
 }
+
+// # of measurement layers hit
+Int_t SolTrack::nmHit()
+{
+        Int_t kmh = 0;
+        Double_t R; Double_t phi; Double_t zz;
+        for (Int_t i = 0; i < fG->Nl(); i++)
+        if (HitLayer(i, R, phi, zz))if (fG->isMeasure(i))kmh++;
+        //
+        return kmh;
+}
+//
+
+
 // List of layers hit with intersections
 Int_t SolTrack::HitList(Int_t *&ihh, Double_t *&rhh, Double_t *&zhh)
@@ -161 +212 @@
   return kmh;
 }
+
+// List of XYZ measurements without any error
+Int_t SolTrack::HitListXYZ(Int_t *&ihh, Double_t *&Xh, Double_t *&Yh, Double_t *&Zh)
+{
+        //
+        // Return lists of hits associated to a track for all measurement layers.
+        // Return value is the total number of measurement hits
+        // kmh = total number of measurement layers hit for given type
+        // ihh = pointer to layer number
+        // Xh, Yh, Zh = X, Y, Z of hit - No measurement error - No multiple scattering
+        //
+        //
+        Int_t kmh = 0;  // Number of measurement layers hit
+        for (Int_t i = 0; i < fG->Nl(); i++)
+        {
+                Double_t R; Double_t phi; Double_t zz;
+                if (HitLayer(i, R, phi, zz)) // Only barrel type layers
+                {
+                        if (fG->isMeasure(i))
+                        {
+                                ihh[kmh] = i;
+                                Xh[kmh] = R*cos(phi);
+                                Yh[kmh] = R*sin(phi);
+                                Zh[kmh] = zz;
+                                kmh++;
+                        }
+                }
+        }
+        //
+        return kmh;
+}
+//
 // Covariance matrix estimation
+//
 void SolTrack::CovCalc(Bool_t Res, Bool_t MS)
 {
-  // Input flags:
-  //        Res = .TRUE. turn on resolution effects/Use standard resolutions
-  //            .FALSE. set all resolutions to 0
-  //        MS  = .TRUE. include Multiple Scattering
-  // Assumptions:
-  // 1. Measurement layers can do one or two measurements
-  // 2. On disks at constant z:
-  //    - Upper side measurement is phi
-  //    - Lower side measurement is R
-
-  // Fill list of layers hit
-  Int_t ntry = 0;
-  Int_t ntrymax = 0;
-  Int_t Nhit = nHit(); // Total number of layers hit
-  Double_t *zhh = new Double_t[Nhit]; // z of hit
-  Double_t *rhh = new Double_t[Nhit]; // r of hit
-  Double_t *dhh = new Double_t[Nhit]; // distance of hit from origin
-  Int_t    *ihh = new Int_t[Nhit];    // true index of layer
-  Int_t kmh; // Number of measurement layers hit
-
-  kmh = HitList(ihh, rhh, zhh); // hit layer list
-  Int_t mTot = 0; // Total number of measurements
-  for (Int_t i = 0; i < Nhit; i++)
-  {
-    dhh[i] = TMath::Sqrt(rhh[i] * rhh[i] + zhh[i] * zhh[i]);
-    if (fG->isMeasure(ihh[i])) mTot += fG->lND(ihh[i]); // Count number of measurements
-  }
-  // Order hit list by increasing distance from origin
-  Int_t *hord = new Int_t[Nhit]; // hit order by increasing distance from origin
-  TMath::Sort(Nhit, dhh, hord, kFALSE); // Order by increasing distance from origin
-  Double_t *zh = new Double_t[Nhit]; // d-ordered z of hit
-  Double_t *rh = new Double_t[Nhit]; // d-ordered r of hit
-  Int_t    *ih = new Int_t[Nhit];    // d-ordered true index of layer
-  for (Int_t i = 0; i < Nhit; i++)
-  {
-    Int_t il = hord[i]; // Hit layer numbering
-    zh[i] = zhh[il];
-    rh[i] = rhh[il];
-    ih[i] = ihh[il];
-  }
-  // Store interdistances and multiple scattering angles
-  Double_t sn2t = 1.0 / (1 + ct()*ct()); //sin^2 theta of track
-  Double_t cs2t = 1.0 - sn2t;            //cos^2 theta
-  Double_t snt = TMath::Sqrt(sn2t);      // sin theta
-  Double_t cst = TMath::Sqrt(cs2t);      // cos theta
-  Double_t px0 = pt() * TMath::Cos(phi0()); // Momentum at minimum approach
-  Double_t py0 = pt() * TMath::Sin(phi0());
-  Double_t pz0 = pt() * ct();
-  TMatrixDSym dik(Nhit);  dik.Zero(); // Distances between layers
-  Double_t *thms = new Double_t[Nhit]; // Scattering angles/plane
-  Double_t *cs = new Double_t[Nhit]; // Cosine of angle with layer normal
-  for (Int_t ii = 0; ii < Nhit; ii++) // Hit layer loop
-  {
-    Int_t i = ih[ii]; // Get true layer number
-    Double_t B = C()*TMath::Sqrt((rh[ii] * rh[ii] - D()*D()) / (1 + 2 * C()*D()));
-    Double_t pxi = px0*(1-2*B*B)-2*py0*B*TMath::Sqrt(1-B*B); // Momentum at scattering layer
-    Double_t pyi = py0*(1-2*B*B)+2*px0*B*TMath::Sqrt(1-B*B);
-    Double_t pzi = pz0;
-    Double_t ArgRp = (rh[ii]*C() + (1 + C() * D())*D() / rh[ii]) / (1 + 2 * C()*D());
-    Double_t phi = phi0() + TMath::ASin(ArgRp);
-    Double_t nx = TMath::Cos(phi); // Barrel layer normal
-    Double_t ny = TMath::Sin(phi);
-    Double_t nz = 0.0;
-    if (fG->lTyp(i) == 2) // this is Z layer
-    {
-      nx = 0.0;
-      ny = 0.0;
-      nz = 1.0;
-    }
-    Double_t corr = (pxi*nx + pyi * ny + pzi * nz) / p();
-    cs[ii] = corr;
-    Double_t Rlf = fG->lTh(i) / (corr*fG->lX0(i)); // Rad. length fraction
-    thms[ii] = 0.0136*TMath::Sqrt(Rlf)*(1.0 + 0.038*TMath::Log(Rlf)) / p(); // MS angle
-    if (!MS)thms[ii] = 0;
-    //
-    for (Int_t kk = 0; kk < ii; kk++) // Fill distances between layers
-    {
-      Double_t Ci = C();
-      dik(ii, kk) = (TMath::ASin(Ci*rh[ii])-TMath::ASin(Ci*rh[kk]))/(Ci*snt);
-      dik(kk, ii) = dik(ii, kk);
-    }
-    // Store momentum components for resolution correction cosines
-    Double_t *pRi = new Double_t[Nhit];
-    pRi[ii] = TMath::Abs(pxi * TMath::Cos(phi) + pyi * TMath::Sin(phi)); // Radial component
-    Double_t *pPhi = new Double_t[Nhit];
-    pPhi[ii] = TMath::Abs(pxi * TMath::Sin(phi) - pyi * TMath::Cos(phi)); // Phi component
-  }
-  // Fill measurement covariance
-  Int_t *mTl = new Int_t[mTot]; // Pointer from measurement number to true layer number
-  TMatrixDSym Sm(mTot); Sm.Zero(); // Measurement covariance
-  TMatrixD Rm(mTot, 5); // Derivative matrix
-  Int_t im = 0;
-  // Fill derivatives and error matrix with MS
-  Double_t AngMax = 90.; Double_t AngMx = AngMax * TMath::Pi() / 180.;
-  Double_t csMin = TMath::Cos(AngMx); // Set maximum angle wrt normal
-  //
-  for (Int_t ii = 0; ii < Nhit; ii++)
-  {
-    Int_t i = ih[ii]; // True layer number
-    Int_t ityp  = fG->lTyp(i); // Layer type Barrel or Z
-    Int_t nmeai = fG->lND(i); // # measurements in layer
-    if (fG->isMeasure(i) && nmeai >0 && cs[ii] > csMin)
-    {
-      Double_t Di    = D(); // Get true track parameters
-      Double_t phi0i = phi0();
-      Double_t Ci    = C();
-      Double_t z0i   = z0();
-      Double_t cti   = ct();
-      //
-      Double_t Ri    = rh[ii];
-      Double_t ArgRp = (Ri*Ci + (1 + Ci * Di)*Di / Ri) / (1 + 2 * Ci*Di);
-      Double_t ArgRz = Ci * TMath::Sqrt((Ri*Ri - Di * Di) / (1 + 2 * Ci*Di));
-      TVectorD dRphi(5); dRphi.Zero(); // R-phi derivatives @ const. R
-      TVectorD dRz(5); dRz.Zero(); // z derivatives @ const. R
-      // Derivative overflow protection
-      Double_t dMin = 0.8;
-      dRphi(0) = (1 - 2 * Ci*Ci*Ri*Ri) /
-        TMath::Max(dMin,TMath::Sqrt(1 - ArgRp * ArgRp)); // D derivative
-      dRphi(1) = Ri; // phi0 derivative
-      dRphi(2) = Ri * Ri /
-        TMath::Max(dMin,TMath::Sqrt(1 - ArgRp * ArgRp)); // C derivative
-      dRphi(3) = 0.0; // z0 derivative
-      dRphi(4) = 0.0; // cot(theta) derivative
-
-      dRz(0) = -cti * Di /
-        (Ri*TMath::Max(dMin,TMath::Sqrt(1 - Ci * Ci*Ri*Ri))); // D
-      dRz(1) = 0.0; // Phi0
-      dRz(2) = cti * (Ri*Ci / TMath::Sqrt(1-Ri*Ri*Ci*Ci) -
-        TMath::ASin(Ri*Ci)) / (Ci*Ci); // C
-      dRz(3) = 1.0; // Z0
-      dRz(4) = TMath::ASin(ArgRz) / Ci; // Cot(theta)
-
-      for (Int_t nmi = 0; nmi < nmeai; nmi++)
-      {
-        mTl[im] = i;
-        Double_t stri = 0;
-        Double_t sig = 0;
-        if (nmi + 1 == 1) // Upper layer measurements
-        {
-          stri = fG->lStU(i); // Stereo angle
-          Double_t csa = TMath::Cos(stri);
-          Double_t ssa = TMath::Sin(stri);
-          sig = fG->lSgU(i); // Resolution
-          if (ityp == 1) // Barrel type layer (Measure R-phi, stereo or z at const. R)
-          {
-            // Almost exact solution (CD<<1, D<<R)
-            Rm(im, 0) = csa * dRphi(0) - ssa * dRz(0); // D derivative
-            Rm(im, 1) = csa * dRphi(1) - ssa * dRz(1); // phi0 derivative
-            Rm(im, 2) = csa * dRphi(2) - ssa * dRz(2); // C derivative
-            Rm(im, 3) = csa * dRphi(3) - ssa * dRz(3); // z0 derivative
-            Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4); // cot(theta) derivative
-          }
-          if (ityp == 2) // Z type layer (Measure phi at const. Z)
-          {
-            Rm(im, 0) = 1.0; // D derivative
-            Rm(im, 1) = rh[ii]; // phi0 derivative
-            Rm(im, 2) = rh[ii] * rh[ii]; // C derivative
-            Rm(im, 3) = 0; // z0 derivative
-            Rm(im, 4) = 0; // cot(theta) derivative
-          }
-        }
-        if (nmi + 1 == 2) // Lower layer measurements
-        {
-          stri = fG->lStL(i); // Stereo angle
-          Double_t csa = TMath::Cos(stri);
-          Double_t ssa = TMath::Sin(stri);
-          sig = fG->lSgL(i); // Resolution
-          if (ityp == 1) // Barrel type layer (measure R-phi, stereo or z at const. R)
-          {
-            // Almost exact solution (CD<<1, D<<R)
-            Rm(im, 0) = csa * dRphi(0) - ssa * dRz(0); // D derivative
-            Rm(im, 1) = csa * dRphi(1) - ssa * dRz(1); // phi0 derivative
-            Rm(im, 2) = csa * dRphi(2) - ssa * dRz(2); // C derivative
-            Rm(im, 3) = csa * dRphi(3) - ssa * dRz(3); // z0 derivative
-            Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4); // cot(theta) derivative
-          }
-          if (ityp == 2) // Z type layer (Measure R at const. z)
-          {
-            Rm(im, 0) = 0; // D derivative
-            Rm(im, 1) = 0; // phi0 derivative
-            Rm(im, 2) = 0; // C derivative
-            Rm(im, 3) = -1.0 / ct(); // z0 derivative
-            Rm(im, 4) = -rh[ii] / ct(); // cot(theta) derivative
-          }
-        }
-        // Derivative calculation completed
-        // Now calculate measurement error matrix
-        Int_t km = 0;
-        for (Int_t kk = 0; kk <= ii; kk++)
-        {
-          Int_t k = ih[kk]; // True layer number
-          Int_t ktyp = fG->lTyp(k); // Layer type Barrel or
-          Int_t nmeak = fG->lND(k); // # measurements in layer
-          if (fG->isMeasure(k) && nmeak > 0 &&cs[kk] > csMin)
-          {
-            for (Int_t nmk = 0; nmk < nmeak; nmk++)
-            {
-              Double_t strk = 0;
-              if (nmk + 1 == 1) strk = fG->lStU(k); // Stereo angle
-              if (nmk + 1 == 2) strk = fG->lStL(k); // Stereo angle
-              if (im == km && Res) Sm(im, km) += sig*sig; // Detector resolution on diagonal
-              //
-              // Loop on all layers below for MS contributions
-              for (Int_t jj = 0; jj < kk; jj++)
-              {
-                Double_t di = dik(ii, jj);
-                Double_t dk = dik(kk, jj);
-                Double_t ms = thms[jj];
-                Double_t msk = ms; Double_t msi = ms;
-                if (ityp == 1) msi = ms / snt; // Barrel
-                else if (ityp == 2) msi = ms / cst; // Disk
-                if (ktyp == 1) msk = ms / snt; // Barrel
-                else if (ktyp == 2) msk = ms / cst; // Disk
-                Double_t ci = TMath::Cos(stri); Double_t si = TMath::Sin(stri);
-                Double_t ck = TMath::Cos(strk); Double_t sk = TMath::Sin(strk);
-                Sm(im, km) += di*dk*(ci*ck*ms*ms + si*sk*msi*msk); // Ms contribution
-              }
-              Sm(km, im) = Sm(im, km);
-              km++;
-            }
-          }
-        }
-        im++; mTot = im;
-      }
-    }
-  }
-  Sm.ResizeTo(mTot, mTot);
-  Rm.ResizeTo(mTot, 5);
-
-  // Calculate covariance from derivatives and measurement error matrix
-  TMatrixDSym DSmInv(mTot); DSmInv.Zero();
-  for (Int_t id = 0; id < mTot; id++) DSmInv(id, id) = 1.0 / TMath::Sqrt(Sm(id, id));
-  TMatrixDSym SmN = Sm.Similarity(DSmInv);  // Normalize diagonal to 1
-  // Protected matrix inversions
-  TDecompChol Chl(SmN);
-  TMatrixDSym SmNinv = SmN;
-  if (Chl.Decompose())
-  {
-    Bool_t OK;
-    SmNinv = Chl.Invert(OK);
-  }
-  else
-  {
-    cout << "SolTrack::CovCalc: Error matrix not positive definite. Recovering ...." << endl;
-    if (ntry < ntrymax)
-    {
-      SmNinv.Print();
-      ntry++;
-    }
-    TMatrixDSym rSmN = MakePosDef(SmN); SmN = rSmN;
-    TDecompChol rChl(SmN);
-    SmNinv = SmN;
-    Bool_t OK = rChl.Decompose();
-    SmNinv    = rChl.Invert(OK);
-  }
-  Sm = SmNinv.Similarity(DSmInv); // Error matrix inverted
-  TMatrixDSym H = Sm.SimilarityT(Rm); // Calculate half Hessian
-  // Normalize before inversion
-  const Int_t Npar = 5;
-  TMatrixDSym DHinv(Npar); DHinv.Zero();
-  for (Int_t i = 0; i < Npar; i++)DHinv(i, i) = 1.0 / TMath::Sqrt(H(i, i));
-  TMatrixDSym Hnrm = H.Similarity(DHinv);
-  // Invert and restore
-  Hnrm.Invert();
-  fCov = Hnrm.Similarity(DHinv);
+        //
+        //
+        // Input flags:
+        //                              Res = .TRUE. turn on resolution effects/Use standard resolutions
+        //                                        .FALSE. set all resolutions to 0
+        //                              MS  = .TRUE. include Multiple Scattering
+        //
+        // Assumptions:
+        // 1. Measurement layers can do one or two measurements
+        // 2. On disks at constant z:
+        //              - Upper side measurement is phi
+        //              - Lower side measurement is R
+        //
+        // Fill list of layers hit
+        //
+        Int_t ntry = 0;
+        Int_t ntrymax = 0;
+        Int_t Nhit = nHit();                                    // Total number of layers hit
+        Double_t *zhh = new Double_t[Nhit];             // z of hit
+        Double_t *rhh = new Double_t[Nhit];             // r of hit
+        Double_t *dhh = new Double_t[Nhit];             // distance of hit from origin
+        Int_t    *ihh = new Int_t[Nhit];                // true index of layer
+        Int_t kmh;                                                              // Number of measurement layers hit
+        //
+        kmh = HitList(ihh, rhh, zhh);                   // hit layer list
+        Int_t mTot = 0;                                                 // Total number of measurements
+        for (Int_t i = 0; i < Nhit; i++)
+        {
+                dhh[i] = TMath::Sqrt(rhh[i] * rhh[i] + zhh[i] * zhh[i]);
+                if (fG->isMeasure(ihh[i])) mTot += fG->lND(ihh[i]);     // Count number of measurements
+        }
+        //
+        // Order hit list by increasing distance from origin
+        //
+        Int_t    *hord = new Int_t[Nhit];               // hit order by increasing distance from origin
+        TMath::Sort(Nhit, dhh, hord, kFALSE);   // Order by increasing distance from origin
+        Double_t *zh = new Double_t[Nhit];              // d-ordered z of hit
+        Double_t *rh = new Double_t[Nhit];              // d-ordered r of hit
+        Int_t    *ih = new Int_t[Nhit];                 // d-ordered true index of layer
+        for (Int_t i = 0; i < Nhit; i++)
+        {
+                Int_t il = hord[i];                                     // Hit layer numbering
+                zh[i] = zhh[il];
+                rh[i] = rhh[il];
+                ih[i] = ihh[il];
+        }
+        //
+        // Store interdistances and multiple scattering angles
+        //
+        Double_t sn2t = 1.0 / (1 + ct()*ct());                  //sin^2 theta of track
+        Double_t cs2t = 1.0 - sn2t;                                             //cos^2 theta
+        Double_t snt = TMath::Sqrt(sn2t);                               // sin theta
+        Double_t cst = TMath::Sqrt(cs2t);                               // cos theta
+        Double_t px0 = pt() * TMath::Cos(phi0());               // Momentum at minimum approach
+        Double_t py0 = pt() * TMath::Sin(phi0());
+        Double_t pz0 = pt() * ct();
+        //
+        TMatrixDSym dik(Nhit);  dik.Zero();             // Distances between layers
+        Double_t *thms = new Double_t[Nhit];    // Scattering angles/plane
+        Double_t *cs = new Double_t[Nhit];              // Cosine of angle with layer normal
+        for (Int_t ii = 0; ii < Nhit; ii++)             // Hit layer loop
+        {
+                Int_t i = ih[ii];                                       // Get true layer number
+                Double_t B = C()*TMath::Sqrt((rh[ii] * rh[ii] - D()*D()) / (1 + 2 * C()*D()));
+                Double_t pxi = px0*(1-2*B*B)-2*py0*B*TMath::Sqrt(1-B*B);                // Momentum at scattering layer
+                Double_t pyi = py0*(1-2*B*B)+2*px0*B*TMath::Sqrt(1-B*B);
+                Double_t pzi = pz0;
+                Double_t ArgRp = (rh[ii]*C() + (1 + C() * D())*D() / rh[ii]) / (1 + 2 * C()*D());
+                Double_t phi = phi0() + TMath::ASin(ArgRp);
+                Double_t nx = TMath::Cos(phi);          // Barrel layer normal
+                Double_t ny = TMath::Sin(phi);
+                Double_t nz = 0.0;
+                if (fG->lTyp(i) == 2)                                                           // this is Z layer
+                {
+                        nx = 0.0;
+                        ny = 0.0;
+                        nz = 1.0;
+                }
+                Double_t corr = (pxi*nx + pyi * ny + pzi * nz) / p();
+                cs[ii] = corr;
+                Double_t Rlf = fG->lTh(i) / (corr*fG->lX0(i));          // Rad. length fraction
+                thms[ii] = 0.0136*TMath::Sqrt(Rlf)*(1.0 + 0.038*TMath::Log(Rlf)) / p();         // MS angle
+                if (!MS)thms[ii] = 0;
+                //
+                for (Int_t kk = 0; kk < ii; kk++)       // Fill distances between layers
+                {
+                        //dik(ii, kk) = TMath::Sqrt(pow(rh[ii] - rh[kk], 2) + pow(zh[ii] - zh[kk], 2));
+                        Double_t Ci = C();
+                        dik(ii, kk) = (TMath::ASin(Ci*rh[ii])-TMath::ASin(Ci*rh[kk]))/(Ci*snt);
+                        dik(kk, ii) = dik(ii, kk);
+                }
+                //
+                // Store momentum components for resolution correction cosines
+                //
+                Double_t *pRi = new Double_t[Nhit];
+                pRi[ii] = TMath::Abs(pxi * TMath::Cos(phi) + pyi * TMath::Sin(phi)); // Radial component
+                Double_t *pPhi = new Double_t[Nhit];
+                pPhi[ii] = TMath::Abs(pxi * TMath::Sin(phi) - pyi * TMath::Cos(phi)); // Phi component
+        }
+        //
+        // Fill measurement covariance
+        //
+        Int_t *mTl = new Int_t[mTot];           // Pointer from measurement number to true layer number
+        TMatrixDSym Sm(mTot); Sm.Zero();        // Measurement covariance
+        TMatrixD Rm(mTot, 5);                           // Derivative matrix
+        Int_t im = 0;
+        //
+        // Fill derivatives and error matrix with MS
+        //
+        Double_t AngMax = 90.; Double_t AngMx = AngMax * TMath::Pi() / 180.;
+        Double_t csMin = TMath::Cos(AngMx);     // Set maximum angle wrt normal
+        //
+        for (Int_t ii = 0; ii < Nhit; ii++)
+        {
+                Int_t i = ih[ii];                                       // True layer number
+                Int_t ityp  = fG->lTyp(i);                      // Layer type Barrel or Z
+                Int_t nmeai = fG->lND(i);                       // # measurements in layer
+                if (fG->isMeasure(i) && nmeai >0 && cs[ii] > csMin)
+                {
+                        Double_t Di    = D();                           // Get true track parameters
+                        Double_t phi0i = phi0();
+                        Double_t Ci    = C();
+                        Double_t z0i   = z0();
+                        Double_t cti   = ct();
+                        //
+                        Double_t Ri    = rh[ii];
+                        Double_t ArgRp = (Ri*Ci + (1 + Ci * Di)*Di / Ri) / (1 + 2 * Ci*Di);
+                        Double_t ArgRz = Ci * TMath::Sqrt((Ri*Ri - Di * Di) / (1 + 2 * Ci*Di));
+                        TVectorD dRphi(5); dRphi.Zero();                // R-phi derivatives @ const. R
+                        TVectorD dRz(5); dRz.Zero();                    // z     derivatives @ const. R
+                        //
+                        // Derivative overflow protection
+                        Double_t dMin = 0.8;
+                        dRphi(0) = (1 - 2 * Ci*Ci*Ri*Ri) /
+                                TMath::Max(dMin,TMath::Sqrt(1 - ArgRp * ArgRp));        // D derivative
+                        dRphi(1) = Ri;                                                                                                          // phi0 derivative
+                        dRphi(2) = Ri * Ri /
+                                TMath::Max(dMin,TMath::Sqrt(1 - ArgRp * ArgRp));                                // C derivative
+                        dRphi(3) = 0.0;                                                                                         // z0 derivative
+                        dRphi(4) = 0.0;                                                                                         // cot(theta) derivative
+                        //
+                        dRz(0) = -cti * Di /
+                                (Ri*TMath::Max(dMin,TMath::Sqrt(1 - Ci * Ci*Ri*Ri)));   // D
+                        dRz(1) = 0.0;                                                                                           // Phi0
+                        dRz(2) = cti * (Ri*Ci / TMath::Sqrt(1-Ri*Ri*Ci*Ci) -            // C
+                                TMath::ASin(Ri*Ci)) / (Ci*Ci);
+                        dRz(3) = 1.0;                                                                                           // Z0
+                        dRz(4) = TMath::ASin(ArgRz) / Ci;                                                       // Cot(theta)
+                        //
+                        for (Int_t nmi = 0; nmi < nmeai; nmi++)
+                        {
+                                mTl[im] = i;
+                                Double_t stri = 0;
+                                Double_t sig = 0;
+                                if (nmi + 1 == 1)               // Upper layer measurements
+                                {
+                                        stri = fG->lStU(i);     // Stereo angle
+                                        Double_t csa = TMath::Cos(stri);
+                                        Double_t ssa = TMath::Sin(stri);
+                                        sig = fG->lSgU(i);      // Resolution
+                                        if (ityp == 1)          // Barrel type layer (Measure R-phi, stereo or z at const. R)
+                                        {
+                                                //
+                                                // Almost exact solution (CD<<1, D<<R)
+                                                Rm(im, 0) = csa * dRphi(0) - ssa * dRz(0);      // D derivative
+                                                Rm(im, 1) = csa * dRphi(1) - ssa * dRz(1);      // phi0 derivative
+                                                Rm(im, 2) = csa * dRphi(2) - ssa * dRz(2);      // C derivative
+                                                Rm(im, 3) = csa * dRphi(3) - ssa * dRz(3);      // z0 derivative
+                                                Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4);      // cot(theta) derivative
+                                        }
+                                        if (ityp == 2)          // Z type layer (Measure phi at const. Z)
+                                        {
+                                                Rm(im, 0) = 1.0;                                        // D derivative
+                                                Rm(im, 1) = rh[ii];                                     // phi0 derivative
+                                                Rm(im, 2) = rh[ii] * rh[ii];            // C derivative
+                                                Rm(im, 3) = 0;                                          // z0 derivative
+                                                Rm(im, 4) = 0;                                          // cot(theta) derivative
+                                        }
+                                }
+                                if (nmi + 1 == 2)                       // Lower layer measurements
+                                {
+                                        stri = fG->lStL(i);             // Stereo angle
+                                        Double_t csa = TMath::Cos(stri);
+                                        Double_t ssa = TMath::Sin(stri);
+                                        sig = fG->lSgL(i);              // Resolution
+                                        if (ityp == 1)                  // Barrel type layer (measure R-phi, stereo or z at const. R)
+                                        {
+                                                //
+                                                // Almost exact solution (CD<<1, D<<R)
+                                                Rm(im, 0) = csa * dRphi(0) - ssa * dRz(0);      // D derivative
+                                                Rm(im, 1) = csa * dRphi(1) - ssa * dRz(1);      // phi0 derivative
+                                                Rm(im, 2) = csa * dRphi(2) - ssa * dRz(2);      // C derivative
+                                                Rm(im, 3) = csa * dRphi(3) - ssa * dRz(3);      // z0 derivative
+                                                Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4);      // cot(theta) derivative
+                                        }
+                                        if (ityp == 2)                  // Z type layer (Measure R at const. z)
+                                        {
+                                                Rm(im, 0) = 0;                                                          // D derivative
+                                                Rm(im, 1) = 0;                                                          // phi0 derivative
+                                                Rm(im, 2) = 0;                                                          // C derivative
+                                                Rm(im, 3) = -1.0 / ct();                                        // z0 derivative
+                                                Rm(im, 4) = -rh[ii] / ct();                                     // cot(theta) derivative
+                                        }
+                                }
+                                // Derivative calculation completed
+                                //
+                                // Now calculate measurement error matrix
+                                //
+                                Int_t km = 0;
+                                for (Int_t kk = 0; kk <= ii; kk++)
+                                {
+                                        Int_t k = ih[kk];                                       // True layer number
+                                        Int_t ktyp = fG->lTyp(k);                       // Layer type Barrel or
+                                        Int_t nmeak = fG->lND(k);                       // # measurements in layer
+                                        if (fG->isMeasure(k) && nmeak > 0 &&cs[kk] > csMin)
+                                        {
+                                                for (Int_t nmk = 0; nmk < nmeak; nmk++)
+                                                {
+                                                        Double_t strk = 0;
+                                                        if (nmk + 1 == 1) strk = fG->lStU(k);   // Stereo angle
+                                                        if (nmk + 1 == 2) strk = fG->lStL(k);   // Stereo angle
+                                                        if (im == km && Res) Sm(im, km) += sig*sig;     // Detector resolution on diagonal
+                                                        //
+                                                        // Loop on all layers below for MS contributions
+                                                        for (Int_t jj = 0; jj < kk; jj++)
+                                                        {
+                                                                Double_t di = dik(ii, jj);
+                                                                Double_t dk = dik(kk, jj);
+                                                                Double_t ms = thms[jj];
+                                                                Double_t msk = ms; Double_t msi = ms;
+                                                                if (ityp == 1) msi = ms / snt;                  // Barrel
+                                                                else if (ityp == 2) msi = ms / cst;             // Disk
+                                                                if (ktyp == 1) msk = ms / snt;                  // Barrel
+                                                                else if (ktyp == 2) msk = ms / cst;             // Disk
+                                                                Double_t ci = TMath::Cos(stri); Double_t si = TMath::Sin(stri);
+                                                                Double_t ck = TMath::Cos(strk); Double_t sk = TMath::Sin(strk);
+                                                                Sm(im, km) += di*dk*(ci*ck*ms*ms + si*sk*msi*msk);                      // Ms contribution
+                                                        }
+                                                        //
+                                                        Sm(km, im) = Sm(im, km);
+                                                        km++;
+                                                }
+                                        }
+                                }
+                                im++; mTot = im;
+                        }
+                }
+        }
+        Sm.ResizeTo(mTot, mTot);
+        Rm.ResizeTo(mTot, 5);
+        //
+        //**********************************************************************
+        // Calculate covariance from derivatives and measurement error matrix  *
+        //**********************************************************************
+        //
+        TMatrixDSym DSmInv(mTot); DSmInv.Zero();
+        for (Int_t id = 0; id < mTot; id++) DSmInv(id, id) = 1.0 / TMath::Sqrt(Sm(id, id));
+        TMatrixDSym SmN = Sm.Similarity(DSmInv);        // Normalize diagonal to 1
+        //SmN.Invert();
+        //
+        // Protected matrix inversions
+        //
+        TDecompChol Chl(SmN);
+        TMatrixDSym SmNinv = SmN;
+        if (Chl.Decompose())
+        {
+                Bool_t OK;
+                SmNinv = Chl.Invert(OK);
+        }
+        else
+        {
+                cout << "SolTrack::CovCalc: Error matrix not positive definite. Recovering ...." << endl;
+                if (ntry < ntrymax)
+                {
+                        SmNinv.Print();
+                        ntry++;
+                }
+                TMatrixDSym rSmN = MakePosDef(SmN); SmN = rSmN;
+                TDecompChol rChl(SmN);
+                SmNinv = SmN;
+                Bool_t OK = rChl.Decompose();
+                SmNinv    = rChl.Invert(OK);
+        }
+        Sm = SmNinv.Similarity(DSmInv);                 // Error matrix inverted
+        TMatrixDSym H = Sm.SimilarityT(Rm);             // Calculate half Hessian
+        // Normalize before inversion
+        const Int_t Npar = 5;
+        TMatrixDSym DHinv(Npar); DHinv.Zero();
+        for (Int_t i = 0; i < Npar; i++)DHinv(i, i) = 1.0 / TMath::Sqrt(H(i, i));
+        TMatrixDSym Hnrm = H.Similarity(DHinv);
+        // Invert and restore
+        Hnrm.Invert();
+        fCov = Hnrm.Similarity(DHinv);
 }
 
@@ -434 +546 @@
 TMatrixDSym SolTrack::MakePosDef(TMatrixDSym NormMat)
 {
-  // Input: symmetric matrix with 1's on diagonal
-  // Output: positive definite matrix with 1's on diagonal
-
-  // Default return value
-  TMatrixDSym rMatN = NormMat;
-  // Check the diagonal
-  Bool_t Check = kFALSE;
-  Int_t Size = NormMat.GetNcols();
-  for (Int_t i = 0; i < Size; i++)if (TMath::Abs(NormMat(i, i) - 1.0)>1.0E-15)Check = kTRUE;
-  if (Check)
-  {
-    cout << "SolTrack::MakePosDef: input matrix doesn't have 1 on diagonal. Abort." << endl;
-    return rMatN;
-  }
-  // Diagonalize matrix
-  TMatrixDSymEigen Eign(NormMat);
-  TMatrixD U = Eign.GetEigenVectors();
-  TVectorD lambda = Eign.GetEigenValues();
-  // Reset negative eigenvalues to small positive value
-  TMatrixDSym D(Size); D.Zero(); Double_t eps = 1.0e-13;
-  for (Int_t i = 0; i < Size; i++)
-  {
-    D(i, i) = lambda(i);
-    if (lambda(i) <= 0) D(i, i) = eps;
-  }
-  // Rebuild matrix
-  TMatrixD Ut(TMatrixD::kTransposed, U);
-  TMatrixD rMat = (U*D)*Ut;       // Now it is positive defite
-  // Restore all ones on diagonal
-  for (Int_t i1 = 0; i1 < Size; i1++)
-  {
-    Double_t rn1 = TMath::Sqrt(rMat(i1, i1));
-    for (Int_t i2 = 0; i2 <= i1; i2++)
-    {
-      Double_t rn2 = TMath::Sqrt(rMat(i2, i2));
-      rMatN(i1, i2) = 0.5*(rMat(i1, i2) + rMat(i2, i1)) / (rn1*rn2);
-      rMatN(i2, i1) = rMatN(i1, i2);
-    }
-  }
-  return rMatN;
-}
+        //
+        // Input: symmetric matrix with 1's on diagonal
+        // Output: positive definite matrix with 1's on diagonal
+        //
+        // Default return value
+        TMatrixDSym rMatN = NormMat;
+        // Check the diagonal
+        Bool_t Check = kFALSE;
+        Int_t Size = NormMat.GetNcols();
+        for (Int_t i = 0; i < Size; i++)if (TMath::Abs(NormMat(i, i) - 1.0)>1.0E-15)Check = kTRUE;
+        if (Check)
+        {
+                cout << "SolTrack::MakePosDef: input matrix doesn ot have 1 on diagonal. Abort." << endl;
+                return rMatN;
+        }
+        //
+        // Diagonalize matrix
+        TMatrixDSymEigen Eign(NormMat);
+        TMatrixD U = Eign.GetEigenVectors();
+        TVectorD lambda = Eign.GetEigenValues();
+        // Reset negative eigenvalues to small positive value
+        TMatrixDSym D(Size); D.Zero(); Double_t eps = 1.0e-13;
+        for (Int_t i = 0; i < Size; i++)
+        {
+                D(i, i) = lambda(i);
+                if (lambda(i) <= 0) D(i, i) = eps;
+        }
+         //Rebuild matrix
+        TMatrixD Ut(TMatrixD::kTransposed, U);
+        TMatrixD rMat = (U*D)*Ut;                               // Now it is positive defite
+        // Restore all ones on diagonal
+        for (Int_t i1 = 0; i1 < Size; i1++)
+        {
+                Double_t rn1 = TMath::Sqrt(rMat(i1, i1));
+                for (Int_t i2 = 0; i2 <= i1; i2++)
+                {
+                        Double_t rn2 = TMath::Sqrt(rMat(i2, i2));
+                        rMatN(i1, i2) = 0.5*(rMat(i1, i2) + rMat(i2, i1)) / (rn1*rn2);
+                        rMatN(i2, i1) = rMatN(i1, i2);
+                }
+        }
+        return rMatN;
+}

external/TrackCovariance/SolTrack.h

@@ -3 +3 @@
 
 #include <TMath.h>
+#include <TVector3.h>
 #include <TMatrixDSym.h>
 
@@ -24 +25 @@
   // Constructors
   SolTrack(Double_t *x, Double_t *p, SolGeom *G);
+        SolTrack(TVector3 x, TVector3 p, SolGeom* G);
   SolTrack(Double_t D, Double_t phi0, Double_t C, Double_t z0, Double_t ct, SolGeom *G);
   // Destructor
@@ -57 +59 @@
   // Track hit management
   Int_t nHit();
+  Int_t nmHit();
   Bool_t HitLayer(Int_t Layer, Double_t &R, Double_t &phi, Double_t &zz);
   Int_t HitList(Int_t *&ihh, Double_t *&rhh, Double_t *&zhh);
+  Int_t HitListXYZ(Int_t *&ihh, Double_t *&Xh, Double_t *&Yh, Double_t *&Zh);
+
   // Make normalized matrix positive definite
   TMatrixDSym MakePosDef(TMatrixDSym NormMat);

modules/TrackCovariance.cc

@@ -21 +21 @@
  *  Smears track parameters according to appropriate covariance matrix.
  *
- *  \authors P. Demin - UCLouvain, Louvain-la-Neuve
+ *  \authors F. Bedeschi - INFN Pisa
+*            P. Demin - UCLouvain, Louvain-la-Neuve
  *           M. Selvaggi - CERN
+ *
  *
  */
@@ -50 +52 @@
 
 TrackCovariance::TrackCovariance() :
-  fGeometry(0), fCovariance(0), fItInputArray(0)
+  fGeometry(0), fCovariance(0), fAcx(0), fItInputArray(0)
 {
   fGeometry = new SolGeom();
@@ -70 +72 @@
   fBz = GetDouble("Bz", 0.0);
   fGeometry->Read(GetString("DetectorGeometry", ""));
+  fNMinHits = GetInt("NMinHits", 6);
 
+  // load geometry
   fCovariance->Calc(fGeometry);
+  fCovariance->SetMinHits(fNMinHits);
+  // load geometry
+  fAcx = fCovariance->AccPnt();
 
   // import input array
-
   fInputArray = ImportArray(GetString("InputArray", "TrackMerger/tracks"));
   fItInputArray = fInputArray->MakeIterator();
@@ -97 +103 @@
   Double_t mass, p, pt, q, ct;
   Double_t dd0, ddz, dphi, dct, dp, dpt, dC;
-  
+
 
   fItInputArray->Reset();
@@ -107 +113 @@
     const TLorentzVector &candidateMomentum = candidate->Momentum;
 
+    if ( !fCovariance->IsAccepted(candidateMomentum.Vect()) ) continue;
+
     mass = candidateMomentum.M();
 
@@ -112 +120 @@
 
     mother    = candidate;
-    candidate = static_cast<Candidate *>(candidate->Clone());
+    candidate = static_cast<Candidate*>(candidate->Clone());
 
     candidate->Momentum.SetVectM(track.GetObsP(), mass);
-    
+
     // converting back to mm
     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;
     candidate->Zd = track.GetObsX().Z()*1e03;
-    
+
     candidate->D0       = track.GetObsPar()[0]*1e03;
     candidate->Phi      = track.GetObsPar()[1];
@@ -142 +150 @@
     candidate->Charge   = q;
 
-    dd0       = TMath::Sqrt(track.GetCov()(0, 0))*1e03; 
-    ddz       = TMath::Sqrt(track.GetCov()(3, 3))*1e03; 
-    dphi      = TMath::Sqrt(track.GetCov()(1, 1)); 
-    dct       = TMath::Sqrt(track.GetCov()(4, 4)); 
+    dd0       = TMath::Sqrt(track.GetCov()(0, 0))*1e03;
+    ddz       = TMath::Sqrt(track.GetCov()(3, 3))*1e03;
+    dphi      = TMath::Sqrt(track.GetCov()(1, 1));
+    dct       = TMath::Sqrt(track.GetCov()(4, 4));
     dpt       = 2 * TMath::Sqrt( track.GetCov()(2, 2))*pt*pt / (0.2998*fBz);
     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;
 
-
     candidate->AddCandidate(mother);
 

modules/TrackCovariance.h

@@ -36 +36 @@
 class SolGeom;
 class SolGridCov;
+class AcceptanceClx;
 
 class TrackCovariance: public DelphesModule
@@ -49 +50 @@
 private:
   Double_t fBz;
+  Int_t fNMinHits;
 
   SolGeom *fGeometry;
   SolGridCov *fCovariance;
+
+  AcceptanceClx *fAcx;
 
   TIterator *fItInputArray; //!