Index: Makefile
===================================================================
--- Makefile	(revision 66b1770d998a4d8c06bb01159db955f51e765bb9)
+++ Makefile	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
@@ -635,4 +635,6 @@
 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,4 +1147,5 @@
 	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) \
Index: cards/delphes_card_IDEA.tcl
===================================================================
--- cards/delphes_card_IDEA.tcl	(revision 66b1770d998a4d8c06bb01159db955f51e765bb9)
+++ cards/delphes_card_IDEA.tcl	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
@@ -1,11 +1,11 @@
 ####################################################################                                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,9 +19,8 @@
   MuonTrackingEfficiency
 
-  ChargedHadronMomentumSmearing
-  ElectronMomentumSmearing
-  MuonMomentumSmearing
-
-  TrackMerger 
+  TrackMergerPre
+  TrackSmearing
+
+  TrackMerger
   Calorimeter
   EFlowMerger
@@ -34,4 +33,5 @@
   ElectronIsolation
 
+  MuonFilter
   MuonEfficiency
   MuonIsolation
@@ -42,5 +42,5 @@
   GenJetFinder
   GenMissingET
-  
+
   FastJetFinder
 
@@ -98,5 +98,5 @@
 }
 
-#	(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,5 +119,5 @@
         (energy < 0.5 && energy >= 0.3) * (abs(eta) <= 3.0)    * (0.65) +
         (energy < 0.3) * (abs(eta) <= 3.0)                     * (0.06)
-    } 
+    }
 }
 
@@ -140,67 +140,218 @@
 }
 
-
-########################################
-# 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,17 +375,17 @@
     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,16 +393,16 @@
     }
 
-    #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,11 +410,11 @@
     }
 
-    # 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,10 +426,10 @@
     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,5 +437,5 @@
     }
 
-    # 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,4 +503,18 @@
 }
 
+
+#################
+# Muon filter
+#################
+
+module PdgCodeFilter MuonFilter {
+  set InputArray Calorimeter/eflowTracks
+  set OutputArray muons
+  set Invert true
+  add PdgCode {13}
+  add PdgCode {-13}
+}
+
+
 #####################
 # Electron efficiency
@@ -363,5 +528,5 @@
 
   # efficiency formula for electrons
-  set EfficiencyFormula {          
+  set EfficiencyFormula {
         (energy < 2.0)                                         * (0.000)+
         (energy >= 2.0) * (abs(eta) <= 0.88)                   * (0.99) +
@@ -388,4 +553,5 @@
 }
 
+
 #################
 # Muon efficiency
@@ -393,5 +559,5 @@
 
 module Efficiency MuonEfficiency {
-  set InputArray MuonMomentumSmearing/muons
+  set InputArray MuonFilter/muons
   set OutputArray muons
 
@@ -399,5 +565,5 @@
 
   # efficiency formula for muons
-  set EfficiencyFormula {                                   
+  set EfficiencyFormula {
         (energy < 2.0)                                         * (0.000)+
         (energy >= 2.0) * (abs(eta) <= 0.88)                   * (0.99) +
@@ -547,5 +713,5 @@
 
   # 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,10 +769,10 @@
 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,17 +782,16 @@
     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
 }
-
Index: external/TrackCovariance/AcceptanceClx.cc
===================================================================
--- external/TrackCovariance/AcceptanceClx.cc	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
+++ external/TrackCovariance/AcceptanceClx.cc	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
@@ -0,0 +1,344 @@
+#include <TVector3.h>
+#include "AcceptanceClx.h"
+#include <iostream>
+//
+// Pt splitting routine
+//
+void AcceptanceClx::VecInsert(Int_t i, Float_t x, TVectorF& Vec)
+{
+	// Insert a new element x in location i+1 of vector Vec
+	//
+	Int_t N = Vec.GetNrows();	// Get vector nitial size
+	Vec.ResizeTo(N + 1);		// Increase size
+	for (Int_t j = N - 1; j > i; j--)Vec(j + 1) = Vec(j);	// Shift all elements above i
+	Vec(i+1) = x;
+}
+//
+void AcceptanceClx::SplitPt(Int_t i, TVectorF &AccPt)
+{
+	Int_t Nrows = fAcc.GetNrows();
+	Int_t Ncols = fAcc.GetNcols();
+	TMatrixF AccMod(Nrows + 1, Ncols); // Size output matrix
+	for (Int_t ith = 0; ith < Ncols; ith++)
+	{
+		AccMod(i + 1, ith) = AccPt(ith);
+		for (Int_t ipt = 0; ipt <= i; ipt++) AccMod(ipt, ith) = fAcc(ipt, ith);
+		for (Int_t ipt = i + 1; ipt < Nrows; ipt++) AccMod(ipt + 1, ith) = fAcc(ipt, ith);
+	}
+	//
+	fAcc.ResizeTo(Nrows + 1, Ncols);
+	fAcc = AccMod;
+}
+//
+// Theta splitting routine
+void AcceptanceClx::SplitTh(Int_t i, TVectorF &AccTh)
+{
+	Int_t Nrows = fAcc.GetNrows();
+	Int_t Ncols = fAcc.GetNcols();
+	TMatrixF AccMod(Nrows, Ncols + 1); // Size output matrix
+	for (Int_t ipt = 0; ipt < Nrows; ipt++)
+	{
+		AccMod(ipt, i + 1) = AccTh(ipt);
+		for (Int_t ith = 0; ith <= i; ith++) AccMod(ipt, ith) = fAcc(ipt, ith);
+		for (Int_t ith = i + 1; ith < Ncols; ith++) AccMod(ipt, ith + 1) = fAcc(ipt, ith);
+	}
+	//
+	fAcc.ResizeTo(Nrows, Ncols + 1);
+	fAcc = AccMod;
+}
+//
+// Constructors
+//
+AcceptanceClx::AcceptanceClx(TString InFile)
+{
+	ReadAcceptance(InFile);
+}
+//
+AcceptanceClx::AcceptanceClx(SolGeom* InGeo)
+{
+	// Initializations
+	//
+	//std::cout << "Entered constructor of AccpeptanceClx" << std::endl;
+	// Setup grid
+	// Start grid parameters
+	//
+	// Pt nodes
+	const Int_t NpPtInp = 5;
+	Float_t PtInit[NpPtInp] = { 0., 1., 10., 100., 250. };
+	TVectorF Pta(NpPtInp, PtInit);
+	Int_t NpPt = Pta.GetNrows();	// Nr. of starting pt points
+	// Theta nodes
+	const Int_t NpThInp = 15;
+	Float_t ThInit[NpThInp] = { 0.,5.,10.,20.,30.,40.,50.,90.,
+	130.,140.,150.,160.,170.,175., 180. };
+	TVectorF Tha(NpThInp, ThInit);
+	Int_t NpTh = Tha.GetNrows();	// Nr. of starting theta points
+	//std::cout << "AcceptanceClv:: Pta and Tha arrays defined" << std::endl;
+	//
+	// Grid splitting parameters
+	Float_t dPtMin = 0.2;		// Grid Pt resolution (GeV)
+	Float_t dThMin = 2.0;		// Grid Theta resolution (degrees)
+	Float_t dAmin  = 1.0;		// Minimum # hits step
+	//
+	fAcc.ResizeTo(NpPt, NpTh);
+	//
+	//
+	// Event loop: fill matrix starting nodes
+	//
+	TVector3 xv(0., 0., 0.);
+	for (Int_t ipt = 0; ipt < NpPt; ipt++)	// Scan pt bins
+	{
+		// Momentum in GeV
+		Double_t pt = Pta(ipt);
+		//
+		for (Int_t ith = 0; ith < NpTh; ith++)	// Scan theta bins
+		{
+			//
+			// Theta in from degrees to radians
+			Double_t th = TMath::Pi() * Tha(ith) / 180.;
+			Double_t pz = pt / TMath::Tan(th);
+			TVector3 tp(pt, 0., pz);
+			//
+			// Get number of measurement hits
+			//
+			SolTrack* gTrk = new SolTrack(xv, tp, InGeo);	// Generated track
+			Int_t Mhits = gTrk->nmHit();			// Nr. Measurement hits
+			fAcc(ipt, ith) = (Float_t)Mhits;
+		}
+	}
+	//
+	// Scan nodes and split if needed
+	//
+	Int_t Nsplits = 1;		// Number of split per iteration
+	Int_t Ncycles = 0;		// Number of iterations
+	Int_t MaxSplits = 200;		// Maximum number of splits
+	Int_t NsplitCnt = 0;
+	//
+	while (Nsplits > 0)
+	{
+		Nsplits = 0;
+		// Scan nodes
+		for (Int_t ipt = 0; ipt < NpPt - 1; ipt++)	// Scan pt bins
+		{
+			for (Int_t ith = 0; ith < NpTh - 1; ith++)	// Scan theta bins
+			{
+				Float_t dAp = TMath::Abs(fAcc(ipt + 1, ith) - fAcc(ipt, ith));
+				Float_t dPt = TMath::Abs(Pta(ipt + 1) - Pta(ipt));
+				//
+				// Pt split
+				if (dPt > dPtMin && dAp > dAmin && Nsplits < MaxSplits) {
+					//std::cout << "Pt(" << ipt << ") = " << Pta(ipt) << ", dAp = " << dAp << std::endl;
+					NsplitCnt++;	// Total splits counter
+					Nsplits++;	// Increase splits/cycle
+					NpPt++;		// Increase #pt points
+					//std::cout << "New pt split: dAp = " << dAp << ", dPt = " << dPt <<
+					//	", Nsplits = " << Nsplits << ", NpPt = " << NpPt << std::endl;
+					Float_t newPt = 0.5 * (Pta(ipt + 1) + Pta(ipt));
+					VecInsert(ipt, newPt, Pta);
+					TVectorF AccPt(NpTh);
+					for (Int_t i = 0; i < NpTh; i++)
+					{
+						Double_t pt = newPt;
+						Double_t th = TMath::Pi() * Tha[i] / 180.;
+						Double_t pz = pt / TMath::Tan(th);
+						TVector3 tp(pt, 0., pz);
+						SolTrack* gTrk = new SolTrack(xv, tp, InGeo);	// Generated track
+						Int_t Mhits = gTrk->nmHit();			// Nr. Measurement hits
+						AccPt(i) = (Float_t)Mhits;
+					}
+					SplitPt(ipt, AccPt);
+					// Completed Pt split
+				}
+				//
+				Float_t dAt = TMath::Abs(fAcc(ipt, ith + 1) - fAcc(ipt, ith));
+				Float_t dTh = TMath::Abs(Tha(ith + 1) - Tha(ith));
+				//
+				// Theta split
+				if (dTh > dThMin && dAt > dAmin && Nsplits < MaxSplits) {
+					//std::cout << "Th(" << ith << ") = " << Tha(ith) << ", dAt = " << dAt << std::endl;
+					NsplitCnt++;	// Total splits counter
+					Nsplits++;	// Increase splits
+					NpTh++;		// Increase #pt points
+					//std::cout << "New th split: dAt = " << dAt << ", dTh = " << dTh <<
+					//	", Nsplits = " << Nsplits << ", NpTh = " << NpTh << std::endl;
+					Float_t newTh = 0.5 * (Tha(ith + 1) + Tha(ith));
+					VecInsert(ith, newTh, Tha);
+					TVectorF AccTh(NpPt);
+					for (Int_t i = 0; i < NpPt; i++)
+					{
+						Double_t pt = Pta(i);
+						Double_t th = TMath::Pi() * newTh / 180.;
+						Double_t pz = pt / TMath::Tan(th);
+						TVector3 tp(pt, 0., pz);
+						SolTrack* gTrk = new SolTrack(xv, tp, InGeo);	// Generated track
+						Int_t Mhits = gTrk->nmHit();			// Nr. Measurement hits
+						AccTh(i) = (Float_t)Mhits;
+					}
+					SplitTh(ith, AccTh);
+					// Theta splits completed
+				}
+			}  // End loop on theta nodes
+		}  // End loop on pt nodes
+		Ncycles++;
+	}			// End loop on iteration
+	//
+	std::cout<<"AcceptanceClx:: Acceptance generation completed after "<<Ncycles<<" cycles"
+		<<" and a total number of "<<NsplitCnt<< " splits"<<std::endl;
+	//
+	// Store final variables and parameters
+	//
+	fNPtNodes = NpPt;
+	Int_t PtCk = Pta.GetNrows();
+	if (PtCk != NpPt)
+		std::cout << "AcceptanceClx:: Error in grid generation NpPt=" << NpPt << ", Pta size= " << PtCk << std::endl;
+	fPtArray.ResizeTo(NpPt);
+	fPtArray = Pta;	// Array of Pt nodes
+	//
+	fNThNodes = NpTh;
+	Int_t ThCk = Tha.GetNrows();
+	if (ThCk != NpTh)
+		std::cout << "AcceptanceClx:: Error in grid generation NpTh=" << NpTh << ", Tha size= " << ThCk << std::endl;
+	fThArray.ResizeTo(NpTh);
+	fThArray = Tha;	// Array of Theta nodes
+	//
+	std::cout << "AcceptanceClx:: Acceptance encoding with " << fNPtNodes
+		<<" pt nodes and "<< fNThNodes <<" theta nodes"<< std::endl;
+	Int_t Nrows = fAcc.GetNrows();
+	Int_t Ncols = fAcc.GetNcols();
+	//std::cout<<"AcceptanceClx:: fAcc rows: "<<Nrows<<", cols: "<<Ncols<<std::endl;
+	//std::cout<<"AcceptanceClx:: Pta size: "<<Pta.GetNrows()<<", Tha size: "<<Tha.GetNrows()<<std::endl;;
+	//std::cout<<"AcceptanceClx:: pt nodes"<<std::endl; fPtArray.Print();
+	//std::cout<<"AcceptanceClx:: th nodes"<<std::endl; fThArray.Print();
+}
+
+// Destructor
+AcceptanceClx::~AcceptanceClx()
+{
+	fNPtNodes = 0;
+	fNThNodes = 0;
+	fAcc.Clear();
+	fPtArray.Clear();
+	fThArray.Clear();
+}
+//
+void AcceptanceClx::WriteAcceptance(TFile *fout)
+{
+	//
+	// Write out data
+	TTree* tree = new TTree("treeAcc", "Acceptance tree");
+	TMatrixF* pntMatF = &fAcc;
+	TVectorF* pntVecP = &fPtArray;
+	TVectorF* pntVecT = &fThArray;
+	tree->Branch("AcceptanceMatrix", "TMatrixF", &pntMatF, 64000, 0);
+	tree->Branch("AcceptancePtVec",  "TVectorF", &pntVecP, 64000, 0);
+	tree->Branch("AcceptanceThVec",  "TVectorF", &pntVecT, 64000, 0);
+	tree->Fill();
+	fout->Write();
+}
+//
+void AcceptanceClx::WriteAcceptance(TString OutFile)
+{
+	//
+	// Write out data
+	TFile* fout = new TFile(OutFile,"RECREATE");
+	WriteAcceptance(fout);
+	fout->Close();
+	delete fout;
+}
+//
+void AcceptanceClx::ReadAcceptance(TString InFile)
+{
+	//
+	// Read in data
+	TFile* f = new TFile(InFile, "READ");
+	//
+	// Import TTree
+	TTree* T = (TTree*)f->Get("treeAcc");
+	//
+	// Get matrix
+	TMatrixF* pAcc = new TMatrixF();
+	T->SetBranchAddress("AcceptanceMatrix", &pAcc);
+	T->GetEntry(0);
+	//
+	fNPtNodes = pAcc->GetNrows();
+	fNThNodes = pAcc->GetNcols();
+	fAcc.ResizeTo(fNPtNodes, fNThNodes);
+	fAcc = *pAcc;
+	//
+	// Get Pt grid
+	TVectorF* pPtArray = new TVectorF();
+	T->SetBranchAddress("AcceptancePtVec", &pPtArray);
+	T->GetEntry(0);
+	fPtArray.ResizeTo(fNPtNodes);
+	fPtArray = *pPtArray;
+	//
+	// Get Theta array
+	TVectorF* pThArray = new TVectorF();
+	T->SetBranchAddress("AcceptanceThVec", &pThArray);
+	T->GetEntry(0);
+	fThArray.ResizeTo(fNThNodes);
+	fThArray = *pThArray;
+	//
+	std::cout << "AcceptanceClx::Read complete: Npt= " << fNPtNodes << ", Nth= " << fNThNodes << std::endl;
+	//
+	f->Close();
+	delete f;
+}
+//
+Double_t AcceptanceClx::HitNumber(Double_t pt, Double_t theta)
+{
+	//
+	// Protect against values out of range
+	//std::cout << "AcceptanceClx::HitNumber: just in pt= " << pt << ", theta = " << theta << std::endl;
+	//std::cout << "AcceptanceClx::HitNumber: ptArr(0)= " << fPtArray(0) << ", thArr(0) = " << fThArray(0) << std::endl;
+	//std::cout << "AcceptanceClx::HitNumber: ptArr(E)= " << fPtArray(fNPtNodes - 1)
+	//	<< ", thArr(E) = " << fThArray(fNThNodes - 1) << std::endl;
+	Float_t eps = 1.0e-4;
+	Float_t pt0 = (Float_t)pt;
+	if (pt0 <= fPtArray(0)) pt0 = fPtArray(0) + eps;
+	else if (pt0 >= fPtArray(fNPtNodes-1)) pt0 = fPtArray(fNPtNodes-1) - eps;
+	Float_t th0 = (Float_t)theta;
+	if (th0 <= fThArray(0)) th0 = fThArray(0) + eps;
+	else if (th0 >= fThArray(fNThNodes - 1)) th0 = fThArray(fNThNodes - 1) - eps;
+	//
+	// Find cell
+	Float_t* parray = fPtArray.GetMatrixArray();
+	Int_t ip = TMath::BinarySearch(fNPtNodes, parray, pt0);
+	Float_t* tarray = fThArray.GetMatrixArray();
+	Int_t it = TMath::BinarySearch(fNThNodes, tarray, th0);
+	//
+	if (ip<0 || ip >fNPtNodes - 2)
+	{
+		std::cout << "Search error: (ip, pt) = (" << ip << ", " << pt << "), pt0 = " << pt0 << std::endl;
+		std::cout << "Search error: pt nodes = " << fNPtNodes
+			<< " , last value = " << fPtArray(fNPtNodes - 1) << std::endl;
+	}
+	if (it<0 || ip >fNThNodes - 2)
+	{
+		std::cout << "Search error: (it, th) = (" << it << ", " << theta << "), th0 = " << th0 << std::endl;
+		std::cout << "Search error: th nodes = " << fNThNodes
+			<< " , last value = " << fThArray(fNThNodes - 1) << std::endl;
+	}
+	//
+	// Bilinear interpolation
+	//
+	Double_t spt = (pt0 - fPtArray(ip)) / (fPtArray(ip + 1) - fPtArray(ip));
+	Double_t sth = (th0 - fThArray(it)) / (fThArray(it + 1) - fThArray(it));
+	Double_t A11 = fAcc(ip, it);
+	Double_t A12 = fAcc(ip, it+1);
+	Double_t A21 = fAcc(ip+1, it);
+	Double_t A22 = fAcc(ip+1, it+1);
+	Double_t A = A11 * (1 - spt) * (1 - sth) + A12 * (1 - spt) * sth +
+		A21 * spt * (1 - sth) + A22 * spt * sth;
+	//
+	return A;
+}
+//
+Double_t AcceptanceClx::HitNum(Double_t *x, Double_t *p) // Theta in degrees
+{
+	Double_t pt = x[0];
+	Double_t th = x[1];
+	//
+	return HitNumber(pt, th);
+}
+//
Index: external/TrackCovariance/AcceptanceClx.h
===================================================================
--- external/TrackCovariance/AcceptanceClx.h	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
+++ external/TrackCovariance/AcceptanceClx.h	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
@@ -0,0 +1,60 @@
+//
+#ifndef G__ACCEPTANCECLX_H
+#define G__ACCEPTANCECLX_H
+//
+#include <TMath.h>
+#include <TVectorF.h>
+#include <TMatrixF.h>
+#include <TString.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TF2.h>
+#include <iostream>
+#include <vector>
+#include "SolGeom.h"
+#include "SolTrack.h"
+//
+// Class to create geometry for solenoid geometry
+
+class AcceptanceClx {
+	//
+	// Class to handle storing and retieving of tracking acceptance
+	//
+private:
+	TMatrixF fAcc;		// Acceptance matrix
+	Int_t fNPtNodes;	// Numer of Pt nodes
+	TVectorF fPtArray;	// Array of Pt nodes
+	Int_t fNThNodes;	// Numer of Theta nodes
+	TVectorF fThArray;	// Array of Theta nodes		(Theta in degrees)
+	//
+	// Service routines
+	void VecInsert(Int_t i, Float_t x, TVectorF& Vec);
+	void SplitPt(Int_t i, TVectorF &AccPt);
+	void SplitTh(Int_t i, TVectorF &AccTh);
+public:
+	//
+	// Constructors
+	AcceptanceClx(SolGeom *InGeo);				// Initialize arrays from geometry
+	AcceptanceClx(TString InFile);				// Initialize from acceptance file
+	// Destructor
+	~AcceptanceClx();
+	//
+	// Accessors
+	TMatrixF* GetAccMatrix() { return &fAcc; }
+	Int_t GetNrPt() { return fNPtNodes; }
+	Int_t GetNrTh() { return fNThNodes; }
+
+	TVectorF* GetPtArray() { return &fPtArray; }
+	TVectorF* GetThArray() { return &fThArray; }
+	//
+	// Read and write
+	void ReadAcceptance(TString InFile);	// Stand alone usage
+	void WriteAcceptance(TString OutFile);	// Stand alone usage
+	void WriteAcceptance(TFile *OutFile);
+	//
+	// Function returning interpolated number of hit measurement layers
+	Double_t HitNumber(Double_t pt, Double_t Theta);	// Theta in degrees
+	Double_t HitNum(Double_t *x, Double_t *p);
+};
+
+#endif
Index: external/TrackCovariance/ObsTrk.cc
===================================================================
--- external/TrackCovariance/ObsTrk.cc	(revision 66b1770d998a4d8c06bb01159db955f51e765bb9)
+++ external/TrackCovariance/ObsTrk.cc	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
@@ -44,4 +44,41 @@
 	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,6 +135,6 @@
 	//
 	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,13 +172,13 @@
 	// 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,5 +216,5 @@
 	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,6 +284,2 @@
 	return cILC;
 }
-
-
-
-	
Index: external/TrackCovariance/ObsTrk.h
===================================================================
--- external/TrackCovariance/ObsTrk.h	(revision 66b1770d998a4d8c06bb01159db955f51e765bb9)
+++ external/TrackCovariance/ObsTrk.h	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
@@ -19,5 +19,5 @@
 	// Prefix Gen marks variables before resolution smearing
 	//
-private:	
+private:
 	Double_t fB;						// Solenoid magnetic field
 	SolGridCov *fGC;					// Covariance matrix grid
@@ -25,5 +25,5 @@
 	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,5 +43,5 @@
 	//
 	TVectorD ParToACTS(TVectorD Par);		// Parameter conversion
-	TMatrixDSym CovToACTS(TMatrixDSym Cov);	// Covariance 
+	TMatrixDSym CovToACTS(TMatrixDSym Cov);	// Covariance
 	//
 	// Conversion to ILC parametrization
@@ -55,5 +55,6 @@
 	// 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();
 	//
Index: external/TrackCovariance/SolGridCov.cc
===================================================================
--- external/TrackCovariance/SolGridCov.cc	(revision 66b1770d998a4d8c06bb01159db955f51e765bb9)
+++ external/TrackCovariance/SolGridCov.cc	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
@@ -3,4 +3,5 @@
 #include <TMath.h>
 #include <TVectorD.h>
+#include <TVector3.h>
 #include <TMatrixDSym.h>
 #include <TDecompChol.h>
@@ -36,4 +37,5 @@
 {
   delete[] fCov;
+  delete fAcc;
 }
 
@@ -58,5 +60,49 @@
     }
   }
-}
+
+// 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)
Index: external/TrackCovariance/SolGridCov.h
===================================================================
--- external/TrackCovariance/SolGridCov.h	(revision 66b1770d998a4d8c06bb01159db955f51e765bb9)
+++ external/TrackCovariance/SolGridCov.h	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
@@ -4,4 +4,5 @@
 #include <TVectorD.h>
 #include <TMatrixDSym.h>
+#include "AcceptanceClx.h"
 
 class SolGeom;
@@ -17,4 +18,6 @@
   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,4 +35,13 @@
   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
+
 };
 
Index: external/TrackCovariance/SolTrack.cc
===================================================================
--- external/TrackCovariance/SolTrack.cc	(revision 66b1770d998a4d8c06bb01159db955f51e765bb9)
+++ external/TrackCovariance/SolTrack.cc	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
@@ -45,4 +45,39 @@
   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,4 +104,6 @@
 SolTrack::~SolTrack()
 {
+	delete[] & fp;
+	delete[] & fpar;
   fCov.Clear();
 }
@@ -132,4 +169,18 @@
   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,272 +212,333 @@
   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,43 +546,45 @@
 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;
+}
Index: external/TrackCovariance/SolTrack.h
===================================================================
--- external/TrackCovariance/SolTrack.h	(revision 66b1770d998a4d8c06bb01159db955f51e765bb9)
+++ external/TrackCovariance/SolTrack.h	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
@@ -3,4 +3,5 @@
 
 #include <TMath.h>
+#include <TVector3.h>
 #include <TMatrixDSym.h>
 
@@ -24,4 +25,5 @@
   // 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,6 +59,9 @@
   // 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);
Index: modules/TrackCovariance.cc
===================================================================
--- modules/TrackCovariance.cc	(revision 66b1770d998a4d8c06bb01159db955f51e765bb9)
+++ modules/TrackCovariance.cc	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
@@ -21,6 +21,8 @@
  *  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,5 +52,5 @@
 
 TrackCovariance::TrackCovariance() :
-  fGeometry(0), fCovariance(0), fItInputArray(0)
+  fGeometry(0), fCovariance(0), fAcx(0), fItInputArray(0)
 {
   fGeometry = new SolGeom();
@@ -70,9 +72,13 @@
   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,5 +103,5 @@
   Double_t mass, p, pt, q, ct;
   Double_t dd0, ddz, dphi, dct, dp, dpt, dC;
-  
+
 
   fItInputArray->Reset();
@@ -107,4 +113,6 @@
     const TLorentzVector &candidateMomentum = candidate->Momentum;
 
+    if ( !fCovariance->IsAccepted(candidateMomentum.Vect()) ) continue;
+
     mass = candidateMomentum.M();
 
@@ -112,8 +120,8 @@
 
     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,5 +138,5 @@
     candidate->Yd = track.GetObsX().Y()*1e03;
     candidate->Zd = track.GetObsX().Z()*1e03;
-    
+
     candidate->D0       = track.GetObsPar()[0]*1e03;
     candidate->Phi      = track.GetObsPar()[1];
@@ -142,8 +150,8 @@
     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,5 +168,4 @@
     candidate->TrackResolution = dp / p;
 
-
     candidate->AddCandidate(mother);
 
Index: modules/TrackCovariance.h
===================================================================
--- modules/TrackCovariance.h	(revision 66b1770d998a4d8c06bb01159db955f51e765bb9)
+++ modules/TrackCovariance.h	(revision 90975be9ee5ec87154c92174683dc36e9878eb5b)
@@ -36,4 +36,5 @@
 class SolGeom;
 class SolGridCov;
+class AcceptanceClx;
 
 class TrackCovariance: public DelphesModule
@@ -49,7 +50,10 @@
 private:
   Double_t fBz;
+  Int_t fNMinHits;
 
   SolGeom *fGeometry;
   SolGridCov *fCovariance;
+
+  AcceptanceClx *fAcx;
 
   TIterator *fItInputArray; //!
