diff --git a/Rich/RichFutureGlobalPID/mathematica/Log(Exp(x)-1).nb b/Rich/RichFutureGlobalPID/mathematica/Log(Exp(x)-1).nb
index c356ae46e2a6e8794ea257eaeed17cdfe9d4d603..6fcd54dce910386bf4dc218b0bf719186ae73be7 100644
--- a/Rich/RichFutureGlobalPID/mathematica/Log(Exp(x)-1).nb
+++ b/Rich/RichFutureGlobalPID/mathematica/Log(Exp(x)-1).nb
@@ -10,10 +10,10 @@
 NotebookFileLineBreakTest
 NotebookFileLineBreakTest
 NotebookDataPosition[       145,          7]
-NotebookDataLength[     91983,       1746]
-NotebookOptionsPosition[     88252,       1652]
-NotebookOutlinePosition[     88791,       1671]
-CellTagsIndexPosition[     88748,       1668]
+NotebookDataLength[     93276,       1764]
+NotebookOptionsPosition[     89546,       1670]
+NotebookOutlinePosition[     90085,       1689]
+CellTagsIndexPosition[     90042,       1686]
 WindowFrame->Normal*)
 
 (* Beginning of Notebook Content *)
@@ -104,8 +104,10 @@ Cell[BoxData["0.001`"], "Output",
    3.7209743612138367`*^9}, 3.72097443610177*^9, 3.72102027097941*^9, 
    3.721020355413443*^9, 3.721020386189633*^9, 3.721020467177422*^9, {
    3.721020541967757*^9, 3.721020567702055*^9}, 3.721036033625578*^9, 
-   3.721036065113412*^9},ExpressionUUID->"8f32baf5-e0bb-4bb1-b9c3-\
-e969c1965264"],
+   3.721036065113412*^9, {3.721108938336089*^9, 3.721108956174171*^9}, 
+   3.721109268145401*^9, {3.721109321627973*^9, 
+   3.721109337418486*^9}},ExpressionUUID->"2ad8cd72-9eda-45c2-955e-\
+80d07e9adc9a"],
 
 Cell[BoxData["0.1`"], "Output",
  CellChangeTimes->{
@@ -141,8 +143,10 @@ Cell[BoxData["0.1`"], "Output",
    3.7209743612138367`*^9}, 3.72097443610177*^9, 3.72102027097941*^9, 
    3.721020355413443*^9, 3.721020386189633*^9, 3.721020467177422*^9, {
    3.721020541967757*^9, 3.721020567702055*^9}, 3.721036033625578*^9, 
-   3.721036065113998*^9},ExpressionUUID->"6e09d45e-0b65-4c2a-8a4f-\
-7062ea710457"],
+   3.721036065113412*^9, {3.721108938336089*^9, 3.721108956174171*^9}, 
+   3.721109268145401*^9, {3.721109321627973*^9, 
+   3.721109337419232*^9}},ExpressionUUID->"46ce6695-8cb9-476f-b78d-\
+061627139e67"],
 
 Cell[BoxData["5.`"], "Output",
  CellChangeTimes->{
@@ -178,8 +182,10 @@ Cell[BoxData["5.`"], "Output",
    3.7209743612138367`*^9}, 3.72097443610177*^9, 3.72102027097941*^9, 
    3.721020355413443*^9, 3.721020386189633*^9, 3.721020467177422*^9, {
    3.721020541967757*^9, 3.721020567702055*^9}, 3.721036033625578*^9, 
-   3.7210360651145563`*^9},ExpressionUUID->"294536fa-79c3-4f40-b8bb-\
-519a0a7e1a9a"]
+   3.721036065113412*^9, {3.721108938336089*^9, 3.721108956174171*^9}, 
+   3.721109268145401*^9, {3.721109321627973*^9, 
+   3.721109337420671*^9}},ExpressionUUID->"acdfa391-78dd-42ca-bb87-\
+ff3e26c2a43f"]
 }, Open  ]],
 
 Cell[CellGroupData[{
@@ -215,8 +221,9 @@ Cell[BoxData[
    3.7209715677794313`*^9, 3.7209715678251543`*^9}, {3.7209721391398993`*^9, 
    3.720972139184471*^9}, {3.720972215542848*^9, 3.7209722450770473`*^9}, {
    3.720972789983202*^9, 3.720972791154159*^9}, {3.720972956060377*^9, 
-   3.7209729803148623`*^9}, {3.720973219409432*^9, 
-   3.720973219488625*^9}},ExpressionUUID->"544dea1a-bfc1-4de2-bf81-\
+   3.7209729803148623`*^9}, {3.720973219409432*^9, 3.720973219488625*^9}, {
+   3.721108929487728*^9, 
+   3.721108953297269*^9}},ExpressionUUID->"544dea1a-bfc1-4de2-bf81-\
 5f4780643639"],
 
 Cell[BoxData[
@@ -286,9 +293,11 @@ Cell[BoxData[
    3.720973674453362*^9}, {3.720974341676323*^9, 3.7209743612920218`*^9}, 
    3.720974436185953*^9, 3.721020271011437*^9, 3.721020355435791*^9, 
    3.7210203862143803`*^9, 3.721020467264583*^9, {3.721020542054976*^9, 
-   3.7210205677782917`*^9}, 3.7210360336460123`*^9, 
-   3.721036065133459*^9},ExpressionUUID->"1ca47f50-4bd6-49cc-b000-\
-307f5926ad16"]
+   3.7210205677782917`*^9}, 3.7210360336460123`*^9, 3.721036065133459*^9, {
+   3.7211089383574266`*^9, 3.721108956200571*^9}, 3.721109268170895*^9, {
+   3.72110932165418*^9, 
+   3.721109337450076*^9}},ExpressionUUID->"d73430d8-f285-4f88-a68a-\
+96d58496759e"]
 }, Open  ]],
 
 Cell[CellGroupData[{
@@ -339,25 +348,12 @@ Cell[BoxData[
    3.720974341743112*^9, 3.720974361358839*^9}, 3.7209744362531548`*^9, 
    3.7210202710288258`*^9, 3.7210203554532413`*^9, 3.721020386231166*^9, 
    3.721020467338324*^9, {3.721020542126025*^9, 3.72102056786279*^9}, 
-   3.721036033663787*^9, 
-   3.721036065152995*^9},ExpressionUUID->"19e2c63b-5b84-47b4-b5da-\
-2106465e625c"]
+   3.721036033663787*^9, 3.721036065152995*^9, {3.721108938385405*^9, 
+   3.721108956220058*^9}, 3.721109268188596*^9, {3.721109321672166*^9, 
+   3.721109337477724*^9}},ExpressionUUID->"ead6837b-2a6a-4dad-a6f5-\
+a9fd6f22c148"]
 }, Open  ]],
 
-Cell[BoxData[
- RowBox[{
-  RowBox[{"fullapprox3", "[", "x_", "]"}], ":=", " ", 
-  RowBox[{
-   RowBox[{"Log", "[", "x", "]"}], " ", "+", " ", 
-   RowBox[{"x", "/", "2"}], " ", "+", " ", 
-   RowBox[{
-    RowBox[{"x", "^", "2"}], "/", "24"}]}]}]], "Input",
- CellChangeTimes->{{3.721020237139719*^9, 3.721020258158375*^9}, {
-  3.721020372671953*^9, 3.721020373247445*^9}, {3.721020455487459*^9, 
-  3.721020463494186*^9}, {3.721036029829385*^9, 
-  3.721036061846005*^9}},ExpressionUUID->"784ab6a2-5c1c-4c0d-9221-\
-8763d738a556"],
-
 Cell[BoxData[
  RowBox[{
   RowBox[{"fullapprox1", "[", "x_", "]"}], ":=", 
@@ -390,6 +386,16 @@ Cell[BoxData[
    3.720973598112768*^9},ExpressionUUID->"5148fb47-f497-40fc-bd3b-\
 df38b48b82d7"],
 
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"fullapprox3", "[", "x_", "]"}], ":=", " ", 
+  RowBox[{
+   RowBox[{"Log", "[", "x", "]"}], " ", "+", " ", 
+   RowBox[{"x", "/", "2"}], " ", "+", " ", 
+   RowBox[{
+    RowBox[{"x", "^", "2"}], "/", "24"}]}]}]], "Input",ExpressionUUID->\
+"35617aa4-1128-4f0e-b95a-6decb2381433"],
+
 Cell[CellGroupData[{
 
 Cell[BoxData[
@@ -472,7 +478,7 @@ QsJE5M2itfMILHh7hfv4BAmB+cmC7XIE/rAoX0qJJ0F1u9HjlLkEWpaTRusT
 SWgJ99NvlyVwNLrBtjyFhBXxwStdaQSaf23z6EsjQTRr9vCiDIE5Nq+E2ukk
 BCzKOdhLJXCo5kP0vkwS3m0QLF7+26aL+1Oyski4k39PIpAm8F89VhK0
        "]]},
-     Annotation[#, "Charting`Private`Tag$11778#1"]& ], 
+     Annotation[#, "Charting`Private`Tag$9374#1"]& ], 
     TagBox[
      {RGBColor[1, 0, 0], AbsoluteThickness[1.6], Opacity[1.], 
       LineBox[CompressedData["
@@ -525,7 +531,7 @@ hcC4xFiRwKIPV1kNCSQkaFqF2soTOGVdvohyioR1seKGFDkCbcrJtaZJJLQV
 mRl2yBI4HtvIKU8hwdeRt8pdhkCroXbvL2kk/Arb/qxQmsA8hzd8vUwSrh+/
 FTwgReBoTX/svmwS5JvTdZf9trnOcEpuLgluuVOfQyQJ/B+eXg0/
        "]]},
-     Annotation[#, "Charting`Private`Tag$11778#2"]& ]}, {}, {}},
+     Annotation[#, "Charting`Private`Tag$9374#2"]& ]}, {}, {}},
   AspectRatio->NCache[GoldenRatio^(-1), 0.6180339887498948],
   Axes->{True, True},
   AxesLabel->{None, None},
@@ -608,9 +614,10 @@ FTwgReBoTX/svmwS5JvTdZf9trnOcEpuLgluuVOfQyQJ/B+eXg0/
    3.720974341959416*^9, 3.720974361456439*^9}, 3.720974436382292*^9, 
    3.72102027119792*^9, 3.721020355506941*^9, 3.721020386289876*^9, 
    3.7210204677186413`*^9, {3.721020542416093*^9, 3.721020568179892*^9}, 
-   3.721036033724605*^9, 
-   3.72103606521496*^9},ExpressionUUID->"5dc7932c-c9cb-4f09-af0b-\
-997d8d1ed23f"]
+   3.721036033724605*^9, 3.72103606521496*^9, {3.721108938502149*^9, 
+   3.7211089562753468`*^9}, 3.7211092683534107`*^9, {3.721109321725006*^9, 
+   3.721109337552956*^9}},ExpressionUUID->"fa0f4d22-ff13-41df-86c5-\
+70d1e2aaa390"]
 }, Open  ]],
 
 Cell[CellGroupData[{
@@ -677,7 +684,7 @@ MYPBu4Kb+t/LvZQbQXoUMQbF3Ecvj8p5fHucZHEJg8gyw1zRrLzHl0Xc8mUM
 tJ6536qYx4DW0J+kvyr/4xbL3kjuBWsuX3q3hoHXNXX9D3Iv0If897R1DDp5
 hNVcuedbY7b71CYGD+3K+Z5SDC4QT69ubcl726pWYfI7/wNv9Ghj
        "]]},
-     Annotation[#, "Charting`Private`Tag$11822#1"]& ], 
+     Annotation[#, "Charting`Private`Tag$9418#1"]& ], 
     TagBox[
      {RGBColor[1, 0, 0], AbsoluteThickness[1.6], Opacity[1.], 
       LineBox[CompressedData["
@@ -721,7 +728,7 @@ GwfNEJBpoXvBcYyAKsN8Qpcv8pt3TajjBJTFNi5O/y3an0/Xej1BwFPdD88c
 gOGEuVnJguiOem6cygYB8lYBQi2R85Nb9/xzk4D0yoihcZG7cpl6Klko+ufU
 8fLHIudWrqd6doeAxeCMeGcBARcp5zZ2dwlw/eOpGyGa8z/UB2GJ
        "]]},
-     Annotation[#, "Charting`Private`Tag$11822#2"]& ]}, {}, {}},
+     Annotation[#, "Charting`Private`Tag$9418#2"]& ]}, {}, {}},
   AspectRatio->NCache[GoldenRatio^(-1), 0.6180339887498948],
   Axes->{True, True},
   AxesLabel->{None, None},
@@ -768,9 +775,10 @@ gOGEuVnJguiOem6cygYB8lYBQi2R85Nb9/xzk4D0yoihcZG7cpl6Klko+ufU
    3.72097434201667*^9, 3.720974361525041*^9}, 3.720974436440476*^9, 
    3.721020271231264*^9, 3.721020355544932*^9, 3.7210203863241177`*^9, 
    3.721020467829204*^9, {3.721020542499977*^9, 3.7210205682909393`*^9}, 
-   3.721036033755262*^9, 
-   3.7210360652474213`*^9},ExpressionUUID->"66276b8d-5ab4-4ba0-968a-\
-1028742a7ce7"]
+   3.721036033755262*^9, 3.7210360652474213`*^9, {3.721108938532302*^9, 
+   3.721108956307886*^9}, 3.721109268378477*^9, {3.72110932176017*^9, 
+   3.721109337594676*^9}},ExpressionUUID->"0f879b9b-663f-4be4-a24d-\
+cc6e7fc9a890"]
 }, Open  ]],
 
 Cell[CellGroupData[{
@@ -850,7 +858,7 @@ ZpsFI0t8L700LmeKQBuzaluyRgjU+UjZOfAngWTCxTKYYwRSDq3pFwgI9KbA
 04kySaDXkxrpxbME4h2LMqxZ8g57Mv2D5wnEznZcFrHk9bU3GpotEkiP9e7L
 sRkCXbPc8EsoXPpxkJdB/CLQ/wBJWmti
        "]]},
-     Annotation[#, "Charting`Private`Tag$11866#1"]& ], 
+     Annotation[#, "Charting`Private`Tag$9462#1"]& ], 
     TagBox[
      {RGBColor[1, 0, 0], AbsoluteThickness[1.6], Opacity[1.], 
       LineBox[CompressedData["
@@ -907,7 +915,7 @@ xDnuvl8inkSmahW3CyYJFK5QtlfgLxL13CSdzk4TqCs1xWImgURkRN3w3NyP
 v4cG/PuTSPRsQim3bIFAo44O1z7dIJGFZ3BA2BKBlAscXnJTSTTUwd2+e/lH
 f9XQMOcWiSKtt8yvrPzwXgGV8tsk+h/C7mQm
        "]]},
-     Annotation[#, "Charting`Private`Tag$11866#2"]& ]}, {}, {}},
+     Annotation[#, "Charting`Private`Tag$9462#2"]& ]}, {}, {}},
   AspectRatio->NCache[GoldenRatio^(-1), 0.6180339887498948],
   Axes->{True, True},
   AxesLabel->{None, None},
@@ -946,8 +954,10 @@ f9XQMOcWiSKtt8yvrPzwXgGV8tsk+h/C7mQm
  CellChangeTimes->{
   3.721020355572196*^9, 3.7210203863617897`*^9, 3.721020467916246*^9, {
    3.7210205425869627`*^9, 3.721020568378347*^9}, 3.721036033791637*^9, 
-   3.721036065280263*^9},ExpressionUUID->"b68d9bbc-88d9-47ff-8214-\
-07ac80f9e003"]
+   3.721036065280263*^9, {3.721108938565785*^9, 3.7211089563455763`*^9}, 
+   3.721109268403881*^9, {3.721109321794969*^9, 
+   3.721109337635182*^9}},ExpressionUUID->"7c266f6f-8edc-4a6e-b76d-\
+8ec85e388ee9"]
 }, Open  ]],
 
 Cell[BoxData[
@@ -1216,7 +1226,7 @@ h6hD6eCU52o6WFfyS8rly+CKsN26+H8ntXt/
        5.1137846302819*^-6}}], 
       LineBox[{{0.0017520034900419211`, 5.1137846302819*^-6}, {
        0.0017525391911699328`, -6.430368902816409*^-6}}]},
-     Annotation[#, "Charting`Private`Tag$11908#1"]& ]}, {}, {}},
+     Annotation[#, "Charting`Private`Tag$9504#1"]& ]}, {}, {}},
   AspectRatio->NCache[GoldenRatio^(-1), 0.6180339887498948],
   Axes->{True, True},
   AxesLabel->{None, None},
@@ -1283,8 +1293,10 @@ h6hD6eCU52o6WFfyS8rly+CKsN26+H8ntXt/
    3.720974361658743*^9}, 3.720974436564753*^9, 3.721020271288372*^9, 
    3.72102035562921*^9, 3.7210203864206057`*^9, 3.721020468252619*^9, {
    3.721020542811964*^9, 3.721020568697794*^9}, 3.7210360338499928`*^9, 
-   3.721036065338706*^9},ExpressionUUID->"21a13aa0-d288-434b-8261-\
-20acabfc347a"]
+   3.721036065338706*^9, {3.7211089386369143`*^9, 3.721108956408988*^9}, 
+   3.7211092684515753`*^9, {3.721109321855029*^9, 
+   3.721109337692778*^9}},ExpressionUUID->"3dce0d1b-d710-49de-8a3e-\
+6103d28092b6"]
 }, Open  ]],
 
 Cell[CellGroupData[{
@@ -1533,7 +1545,7 @@ O6Jy4f8AXMpMqA==
        0.5223833909925699, -0.00004533017015007079}, {
        0.5256190478570228, -0.00005990932286953354}, {
        0.5286746903331817, -0.00007305880620850901}}]},
-     Annotation[#, "Charting`Private`Tag$11949#1"]& ]}, {}, {}},
+     Annotation[#, "Charting`Private`Tag$9545#1"]& ]}, {}, {}},
   AspectRatio->NCache[GoldenRatio^(-1), 0.6180339887498948],
   Axes->{True, True},
   AxesLabel->{None, None},
@@ -1580,8 +1592,10 @@ O6Jy4f8AXMpMqA==
    3.720974361732049*^9}, 3.7209744366207857`*^9, 3.721020271328685*^9, 
    3.7210203556640882`*^9, 3.721020386454915*^9, 3.721020468340653*^9, {
    3.7210205429025917`*^9, 3.721020568808782*^9}, 3.72103603388811*^9, 
-   3.721036065373519*^9},ExpressionUUID->"207ed403-2b5b-4515-be4d-\
-572a16114503"]
+   3.721036065373519*^9, {3.72110893867492*^9, 3.721108956450013*^9}, 
+   3.721109268488559*^9, {3.721109321892415*^9, 
+   3.721109337730337*^9}},ExpressionUUID->"332b2707-35a2-4713-ae64-\
+7a8aa8921972"]
 }, Open  ]],
 
 Cell[CellGroupData[{
@@ -1619,8 +1633,10 @@ Cell["\<\
    3.720974342308324*^9, 3.7209743617917557`*^9}, 3.720974436666292*^9, 
    3.72102027133853*^9, 3.721020355673856*^9, 3.7210203864643707`*^9, 
    3.721020468415882*^9, {3.721020542998926*^9, 3.7210205688807898`*^9}, 
-   3.721036065383286*^9},ExpressionUUID->"d4dbbe45-fba7-49fb-82e3-\
-1c198fc299b4"]
+   3.721036065383286*^9, {3.72110893868503*^9, 3.7211089564599657`*^9}, 
+   3.721109268498831*^9, {3.721109321902131*^9, 
+   3.7211093377401323`*^9}},ExpressionUUID->"f587a756-a889-4225-8fa9-\
+e74b23fa4446"]
 }, Open  ]],
 
 Cell[CellGroupData[{
@@ -1645,13 +1661,15 @@ Cell["\<\
    3.720973631646098*^9, 3.7209736751325397`*^9}, {3.720974342372404*^9, 
    3.720974361861771*^9}, 3.720974436732634*^9, 3.721020271356639*^9, 
    3.721020355689389*^9, 3.721020386483366*^9, 3.721020468503758*^9, {
-   3.721020543070097*^9, 3.721020568963121*^9}, 
-   3.7210360654000587`*^9},ExpressionUUID->"1792d1c8-ce5a-4972-b737-\
-966ffdfba480"]
+   3.721020543070097*^9, 3.721020568963121*^9}, 3.7210360654000587`*^9, {
+   3.72110893870453*^9, 3.7211089564783163`*^9}, 3.721109268518867*^9, {
+   3.721109321917596*^9, 
+   3.7211093377550383`*^9}},ExpressionUUID->"3c20e4ab-b26e-4d61-be59-\
+9392931f1821"]
 }, Open  ]]
 },
 WindowSize->{971, 1002},
-WindowMargins->{{316, Automatic}, {4, Automatic}},
+WindowMargins->{{320, Automatic}, {4, Automatic}},
 FrontEndVersion->"11.1 for Linux x86 (64-bit) (April 18, 2017)",
 StyleDefinitions->Notebook[{
    Cell[
@@ -1677,76 +1695,76 @@ Cell[1416, 38, 280, 7, 55, "Input", "ExpressionUUID" -> \
 Cell[CellGroupData[{
 Cell[1721, 49, 1394, 21, 77, "Input", "ExpressionUUID" -> \
 "70764c1e-89d9-4fbb-8b21-680ebdd742fc"],
-Cell[3118, 72, 2476, 35, 32, "Output", "ExpressionUUID" -> \
-"8f32baf5-e0bb-4bb1-b9c3-e969c1965264"],
-Cell[5597, 109, 2474, 35, 32, "Output", "ExpressionUUID" -> \
-"6e09d45e-0b65-4c2a-8a4f-7062ea710457"],
-Cell[8074, 146, 2475, 35, 32, "Output", "ExpressionUUID" -> \
-"294536fa-79c3-4f40-b8bb-519a0a7e1a9a"]
+Cell[3118, 72, 2598, 37, 32, "Output", "ExpressionUUID" -> \
+"2ad8cd72-9eda-45c2-955e-80d07e9adc9a"],
+Cell[5719, 111, 2596, 37, 32, "Output", "ExpressionUUID" -> \
+"46ce6695-8cb9-476f-b78d-061627139e67"],
+Cell[8318, 150, 2595, 37, 32, "Output", "ExpressionUUID" -> \
+"acdfa391-78dd-42ca-bb87-ff3e26c2a43f"]
 }, Open  ]],
 Cell[CellGroupData[{
-Cell[10586, 186, 2063, 33, 34, "Input", "ExpressionUUID" -> \
+Cell[10950, 192, 2113, 34, 34, "Input", "ExpressionUUID" -> \
 "544dea1a-bfc1-4de2-bf81-5f4780643639"],
-Cell[12652, 221, 4284, 69, 55, "Output", "ExpressionUUID" -> \
-"1ca47f50-4bd6-49cc-b000-307f5926ad16"]
+Cell[13066, 228, 4407, 71, 55, "Output", "ExpressionUUID" -> \
+"d73430d8-f285-4f88-a68a-96d58496759e"]
 }, Open  ]],
 Cell[CellGroupData[{
-Cell[16973, 295, 617, 14, 34, "Input", "ExpressionUUID" -> \
+Cell[17510, 304, 617, 14, 34, "Input", "ExpressionUUID" -> \
 "fa1218a9-ca1e-4ca1-8f39-be7d54d3acab"],
-Cell[17593, 311, 1615, 32, 55, "Output", "ExpressionUUID" -> \
-"19e2c63b-5b84-47b4-b5da-2106465e625c"]
+Cell[18130, 320, 1733, 33, 55, "Output", "ExpressionUUID" -> \
+"ead6837b-2a6a-4dad-a6f5-a9fd6f22c148"]
 }, Open  ]],
-Cell[19223, 346, 522, 12, 34, "Input", "ExpressionUUID" -> \
-"784ab6a2-5c1c-4c0d-9221-8763d738a556"],
-Cell[19748, 360, 872, 17, 34, "Input", "ExpressionUUID" -> \
+Cell[19878, 356, 872, 17, 34, "Input", "ExpressionUUID" -> \
 "0f0e30ff-5c8e-4d65-a7d7-1494d34243e5"],
-Cell[20623, 379, 448, 11, 34, "Input", "ExpressionUUID" -> \
+Cell[20753, 375, 448, 11, 34, "Input", "ExpressionUUID" -> \
 "5148fb47-f497-40fc-bd3b-df38b48b82d7"],
+Cell[21204, 388, 309, 8, 34, "Input", "ExpressionUUID" -> \
+"35617aa4-1128-4f0e-b95a-6decb2381433"],
 Cell[CellGroupData[{
-Cell[21096, 394, 1241, 24, 34, "Input", "ExpressionUUID" -> \
+Cell[21538, 400, 1241, 24, 34, "Input", "ExpressionUUID" -> \
 "f4a96435-b8fb-4c6c-8b82-8a9520e6e33a"],
-Cell[22340, 420, 11006, 192, 290, "Output", "ExpressionUUID" -> \
-"5dc7932c-c9cb-4f09-af0b-997d8d1ed23f"]
+Cell[22782, 426, 11126, 193, 290, "Output", "ExpressionUUID" -> \
+"fa0f4d22-ff13-41df-86c5-70d1e2aaa390"]
 }, Open  ]],
 Cell[CellGroupData[{
-Cell[33383, 617, 595, 15, 34, "Input", "ExpressionUUID" -> \
+Cell[33945, 624, 595, 15, 34, "Input", "ExpressionUUID" -> \
 "82366ad3-f2a6-4cfa-b93d-d710ea2a381f"],
-Cell[33981, 634, 7258, 138, 284, "Output", "ExpressionUUID" -> \
-"66276b8d-5ab4-4ba0-968a-1028742a7ce7"]
+Cell[34543, 641, 7373, 139, 284, "Output", "ExpressionUUID" -> \
+"0f879b9b-663f-4be4-a24d-cc6e7fc9a890"]
 }, Open  ]],
 Cell[CellGroupData[{
-Cell[41276, 777, 547, 15, 34, "Input", "ExpressionUUID" -> \
+Cell[41953, 785, 547, 15, 34, "Input", "ExpressionUUID" -> \
 "019efc88-b53f-48a0-83c2-b157e309b7c5"],
-Cell[41826, 794, 8134, 155, 300, "Output", "ExpressionUUID" -> \
-"b68d9bbc-88d9-47ff-8214-07ac80f9e003"]
+Cell[42503, 802, 8256, 157, 300, "Output", "ExpressionUUID" -> \
+"7c266f6f-8edc-4a6e-b76d-8ec85e388ee9"]
 }, Open  ]],
-Cell[49975, 952, 624, 13, 34, "Input", "ExpressionUUID" -> \
+Cell[50774, 962, 624, 13, 34, "Input", "ExpressionUUID" -> \
 "7e13f64f-1c1d-4d40-8972-edd78cf14e7c"],
-Cell[50602, 967, 330, 9, 34, "Input", "ExpressionUUID" -> \
+Cell[51401, 977, 330, 9, 34, "Input", "ExpressionUUID" -> \
 "62a0911c-bdac-47d1-8aa9-2f66e7220c55"],
 Cell[CellGroupData[{
-Cell[50957, 980, 1068, 18, 34, "Input", "ExpressionUUID" -> \
+Cell[51756, 990, 1068, 18, 34, "Input", "ExpressionUUID" -> \
 "bbcc344a-5e17-4d62-9e00-32e28e82d17d"],
-Cell[52028, 1000, 16607, 286, 257, "Output", "ExpressionUUID" -> \
-"21a13aa0-d288-434b-8261-20acabfc347a"]
+Cell[52827, 1010, 16732, 288, 257, "Output", "ExpressionUUID" -> \
+"3dce0d1b-d710-49de-8a3e-6103d28092b6"]
 }, Open  ]],
 Cell[CellGroupData[{
-Cell[68672, 1291, 302, 8, 34, "Input", "ExpressionUUID" -> \
+Cell[69596, 1303, 302, 8, 34, "Input", "ExpressionUUID" -> \
 "4198608e-3e77-4484-bcce-7b1cc34384ea"],
-Cell[68977, 1301, 16091, 282, 265, "Output", "ExpressionUUID" -> \
-"207ed403-2b5b-4515-be4d-572a16114503"]
+Cell[69901, 1313, 16211, 284, 265, "Output", "ExpressionUUID" -> \
+"332b2707-35a2-4713-ae64-7a8aa8921972"]
 }, Open  ]],
 Cell[CellGroupData[{
-Cell[85105, 1588, 190, 4, 32, "Input", "ExpressionUUID" -> \
+Cell[86149, 1602, 190, 4, 32, "Input", "ExpressionUUID" -> \
 "653eef77-941c-4885-994b-11d2f55f5324"],
-Cell[85298, 1594, 1853, 28, 103, "Output", "ExpressionUUID" -> \
-"d4dbbe45-fba7-49fb-82e3-1c198fc299b4"]
+Cell[86342, 1608, 1978, 30, 103, "Output", "ExpressionUUID" -> \
+"f587a756-a889-4225-8fa9-e74b23fa4446"]
 }, Open  ]],
 Cell[CellGroupData[{
-Cell[87188, 1627, 192, 4, 32, "Input", "ExpressionUUID" -> \
+Cell[88357, 1643, 192, 4, 32, "Input", "ExpressionUUID" -> \
 "f34c1945-b2e5-4254-98a2-85e29a2a395a"],
-Cell[87383, 1633, 853, 16, 84, "Output", "ExpressionUUID" -> \
-"1792d1c8-ce5a-4972-b737-966ffdfba480"]
+Cell[88552, 1649, 978, 18, 84, "Output", "ExpressionUUID" -> \
+"3c20e4ab-b26e-4d61-be59-9392931f1821"]
 }, Open  ]]
 }
 ]
diff --git a/Rich/RichFutureGlobalPID/scripts/drawfuncs.C b/Rich/RichFutureGlobalPID/scripts/drawfuncs.C
new file mode 100644
index 0000000000000000000000000000000000000000..dbe89ab38a0d29e111093fa209c8e9d6812b4d3b
--- /dev/null
+++ b/Rich/RichFutureGlobalPID/scripts/drawfuncs.C
@@ -0,0 +1,88 @@
+
+void drawfuncs()
+{
+
+  auto c = new TCanvas("log(exp(x)-1)","c",1000,800);
+
+  const double fmin = 0.001;
+  const double fmax = 5.0;
+
+  // the baseline using STD functions
+  auto stdf = []( double*x, double * ) 
+    {
+      return std::log( std::exp( x[0] ) - 1.0 ); 
+    };
+  auto tf1 = new TF1( "log(exp(x)-1)", stdf, fmin, fmax, 1 );
+  tf1->SetLineColor(kBlack);
+  tf1->Draw();
+
+  // approximation
+  auto a1 = []( double * xp, double * ) 
+    {
+      // using log( e^x - 1 ) ~= log(x) + x/2 + x^2/24 - x^4/2880 + x^6/181440
+
+      const float x = (float)xp[0];
+
+      // The log part
+      const union { float f; std::uint32_t i; } vx = { x };
+      const union { std::uint32_t i; float f; } mx = { ( vx.i & 0x007FFFFF ) | 0x3f000000 };
+      const auto y = float(vx.i) * 1.1920928955078125e-7f;
+      auto z = 0.69314718f * ( y - 124.22551499f - 1.498030302f * mx.f 
+                               - 1.72587999f / ( 0.3520887068f + mx.f ) );
+
+      // the power series
+      const float a( 1.0/24.0 );
+      const float b( 0.5      );
+      z += ( ( ( a * x ) + b ) * x ); // x and x^2 terms
+      z -= std::pow(x,4)/2800.0;      // x^4 term
+      z += std::pow(x,6)/181440.0;    // x^6 term
+      
+      // return
+      return z;
+    };
+  auto tf2 = new TF1( "approx1", a1, fmin, fmax, 1 );
+  tf2->SetLineColor(kRed);
+  tf2->Draw("SAME");
+
+  // faster approximation
+  auto a2 = []( double * xp, double * ) 
+    {
+      // using log( e^x - 1 ) ~= log(x) + x/2 + x^2/24 - x^4/2880 + x^6/181440
+
+      const float x = (float)xp[0];
+
+      // The log part
+      const union { float f; std::uint32_t i; } vx = { x };
+      auto z = ( float(vx.i) * 8.2629582881927490e-8f ) - 87.989971088f;
+
+      // the power series
+      const float a( 1.0/24.0 );
+      const float b( 0.5      );
+      z += b*x;   // x term
+      z += a*x*x; // x^2 term 
+      
+      // return
+      return z;
+    };
+  auto tf3 = new TF1( "approx2", a2, fmin, fmax, 1 );
+  tf3->SetLineColor(kBlue);
+  tf3->Draw("SAME");
+
+  c = new TCanvas("log(exp(x)-1) log10(diff)","c",1000,800);
+
+  auto diff1 = [&stdf,&a1]( double * x, double * p ) { return std::log10( std::fabs( a1(x,p) - stdf(x,p) ) ); };
+  auto diff2 = [&stdf,&a2]( double * x, double * p ) { return std::log10( std::fabs( a2(x,p) - stdf(x,p) ) ); };
+
+  auto tdiff1 = new TF1( "log(exp(x)-1) | log10( fabs(approx-STD) )", diff1, fmin, fmax, 1 );
+  auto tdiff2 = new TF1( "diff2",                                     diff2, fmin, fmax, 1 );
+
+  tdiff1->SetLineColor(kRed);
+  tdiff2->SetLineColor(kBlue);
+
+  tdiff1->SetMinimum(-8);
+  tdiff1->SetMaximum(0);
+
+  tdiff1->Draw();
+  tdiff2->Draw("SAME");
+
+}
diff --git a/Rich/RichFutureGlobalPID/src/RichGlobalPIDLikelihoodMinimiser.cpp b/Rich/RichFutureGlobalPID/src/RichGlobalPIDLikelihoodMinimiser.cpp
index 3d5e2ea496b0b49f226d12d757ef64aa8656eb0b..49f6462521cb154c7da8c764db66e4566f14ce54 100644
--- a/Rich/RichFutureGlobalPID/src/RichGlobalPIDLikelihoodMinimiser.cpp
+++ b/Rich/RichFutureGlobalPID/src/RichGlobalPIDLikelihoodMinimiser.cpp
@@ -33,9 +33,6 @@ StatusCode LikelihoodMinimiser::initialize()
   auto sc = MultiTransformer::initialize();
   if ( !sc ) return sc;
 
-  // Initialise the interpolator
-  m_logExpLookUp.init( 1.0 ); // the max value
-
   // Initialise parameters
   m_logMinSig = logExp(m_minSig);
 
diff --git a/Rich/RichFutureGlobalPID/src/RichGlobalPIDLikelihoodMinimiser.h b/Rich/RichFutureGlobalPID/src/RichGlobalPIDLikelihoodMinimiser.h
index 24db1009a3e059b1b8d5a6f44624a3f3cf46d40f..32ed22f75821761b9bf9141e5ec1507d50f56522 100644
--- a/Rich/RichFutureGlobalPID/src/RichGlobalPIDLikelihoodMinimiser.h
+++ b/Rich/RichFutureGlobalPID/src/RichGlobalPIDLikelihoodMinimiser.h
@@ -303,7 +303,11 @@ namespace Rich
           /// Prefered implementation of log( e^x -1 )
           inline FloatType logExp( const FloatType x ) const noexcept
           {
-            return m_logExpLookUp.lookup(x);
+            using namespace Rich::Maths::Approx;
+            const FloatType a( 1.0/24.0 );
+            const FloatType b( 0.5      );
+            //return approx_log(x) + ( ( ( a * x ) + b ) * x );
+            return vapprox_log(x) + ( ( ( a * x ) + b ) * x );
           }
 
           /// log( exp(x) - 1 ) or an approximation for small signals
@@ -399,9 +403,6 @@ namespace Rich
 
           /// Cached value of log( exp(m_minSig) - 1 ) for efficiency
           FloatType m_logMinSig{0};
-          
-          /// Look up table for log(exp(x)-1) function
-          LogExpLookUp<FloatType,10000> m_logExpLookUp;
 
           /// Minimum signal value for full calculation of log(exp(signal)-1)
           Gaudi::Property<FloatType> m_minSig
diff --git a/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.cpp b/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.cpp
index 54e089b3db08577de6f00c748885971e2d08a1cb..85bcb1ad71bc64ec169bac16ff8ccd2d413c8f8c 100644
--- a/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.cpp
+++ b/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.cpp
@@ -35,11 +35,11 @@ StatusCode SIMDLikelihoodMinimiser::initialize()
   if ( !sc ) return sc;
 
   // Initialise the interpolator for f = log ( exp(x) - 1 )
-  m_logExpLookUp.init( FloatType(5), 
-                       []( const FloatType x )
-                       { const FloatType a = vdt::fast_exp(x) - 1.0;
-                         const FloatType b = std::max(a,std::numeric_limits<FloatType>::min());
-                         return vdt::fast_log(b); } );
+  //m_logExpLookUp.init( FloatType(5), 
+  //                     []( const FloatType x )
+  //                     { const FloatType a = vdt::fast_exp(x) - 1.0;
+  //                       const FloatType b = std::max(a,std::numeric_limits<FloatType>::min());
+  //                       return vdt::fast_log(b); } );
                          
   // Initialise parameters
   m_minSigSIMD = SIMDFP( m_minSig );
diff --git a/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.h b/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.h
index 03708a0a5f94d0afbf00c34568b0de66cf1b2e89..868f7092183d2c9c8074ff930209736762398bee 100644
--- a/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.h
+++ b/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.h
@@ -35,7 +35,7 @@
 #include "RichUtils/ZipRange.h"
 #include "RichUtils/RichMap.h"
 #include "RichUtils/FastMaths.h"
-#include "RichUtils/LookupTableInterpolator.h"
+//#include "RichUtils/LookupTableInterpolator.h"
 
 namespace Rich 
 {
@@ -194,7 +194,83 @@ namespace Rich
           template< typename TYPE >
           inline TYPE full_logExp( const TYPE & x ) const noexcept
           {
-            return vdt::fast_log( vdt::fast_exp(x) - 1.0 );
+            using namespace Rich::Maths;
+            return fast_log( fast_exp(x) - TYPE(1.0) );
+          }
+
+          /// Use rational interpolation
+          /// Use two functions, one for 0.001 to 0.1 and a second for 0.1 to 5
+          template< typename TYPE,
+                    typename std::enable_if<!std::is_arithmetic<TYPE>::value>::type * = nullptr >
+          inline TYPE rational_interp_logExp( const TYPE & x ) const noexcept
+          {
+
+            // test which interpolators are required
+            const auto m = ( x > TYPE(0.1) );
+
+            // the return value
+            TYPE res;
+
+            // 0.001 to 0.1
+            // top params
+            if ( !all_of(m) )
+            {
+              const TYPE t4( 1.5368801356135493e6 );
+              const TYPE t3( -923340.7133801008   );
+              const TYPE t2( -180149.47011640694  );
+              const TYPE t1( -3763.9820937450886  );
+              const TYPE t0( -8.153881974971727   );
+              // bottom params
+              const TYPE b4( 765067.9490218208  );
+              const TYPE b3( 631524.6174061331  );
+              const TYPE b2( 53227.905690870386 );
+              const TYPE b1( 710.6892499540406  );
+              const TYPE b0( 1 );
+              // compute ratio ( FMA friendly )
+              const auto t = ( ( ( ( ( ( ( t4 * x ) + t3 ) * x ) + t2 ) * x ) + t1 ) * x ) + t0;
+              const auto b = ( ( ( ( ( ( ( b4 * x ) + b3 ) * x ) + b2 ) * x ) + b1 ) * x ) + b0;
+              // set res
+              res = t/b;
+            }
+           
+            // 0.1 to 5.0
+            if ( any_of(m) )
+            {
+              // top params
+              const TYPE n3( 17.67838014763876   );
+              const TYPE n2( 34.24320868914387   );
+              const TYPE n1( -25.594854219076243 );
+              const TYPE n0( -4.607325499367173  );
+              // bottom parms
+              const TYPE d4( 0.019476828793646282 );
+              const TYPE d3( -0.4365084410635895  );
+              const TYPE d2( 21.554896658342326   );
+              const TYPE d1( 17.994362235519787   );
+              const TYPE d0( 1 );
+              // compute ratio ( FMA friendly )
+              const auto n = ( ( ( ( ( (              n3 ) * x ) + n2 ) * x ) + n1 ) * x ) + n0;
+              const auto d = ( ( ( ( ( ( ( d4 * x ) + d3 ) * x ) + d2 ) * x ) + d1 ) * x ) + d0;
+              // set res for those required
+              res(m) = n/d;
+            }
+  
+            // return
+            return res;
+          }     
+          
+          /// Approximate log( e^x - 1 ) for SIMD types
+          template< typename TYPE,
+                    typename std::enable_if<!std::is_arithmetic<TYPE>::value>::type * = nullptr >
+          inline TYPE approx_logExp( const TYPE & x ) const noexcept
+          {
+            using namespace Rich::Maths::Approx;
+            // Use power series expansion
+            // log( e^x - 1 ) ~= log(x) + x/2 + x^2/24
+            // works well for x ~ 0.001 to 5 
+            const TYPE a( 1.0/24.0 );
+            const TYPE b( 0.5      );
+            //return approx_log(x) + ( ( ( a * x ) + b ) * x );
+            return vapprox_log(x) + ( ( ( a * x ) + b ) * x );
           }
 
           /// Fast implementation of log( e^x - 1 ) for SIMD types
@@ -203,10 +279,16 @@ namespace Rich
           inline TYPE fast_logExp( const TYPE & x ) const noexcept
           {
             // Use the fast VDT inspired methods
-            //using namespace Rich::SIMD::Maths;
-            //return fast_log( fast_exp(x) - TYPE::One() );
+            //return full_logExp(x);
+
+            // Use rational interpolation
+            //return rational_interp_logExp(x);
+
+            // Use the fast approoximation
+            return approx_logExp(x);
+
             // use interpolator
-            return m_logExpLookUp(x);
+            //return m_logExpLookUp(x);
           }
 
           /// log( exp(x) - 1 ) or an approximation for small signals
@@ -314,7 +396,7 @@ namespace Rich
           SIMDFP m_logMinSig = SIMDFP::Zero();
           
           /// Look up table for log(exp(x)-1) function
-          LookupTableInterpolatorFromZero<FloatType,10000> m_logExpLookUp;
+          //LookupTableInterpolatorFromZero<FloatType,10000> m_logExpLookUp;
 
           /// Minimum signal value for full calculation of log(exp(signal)-1)
           Gaudi::Property<FloatType> m_minSig
diff --git a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecCherenkovAngles.h b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecCherenkovAngles.h
index 1447eba3a28aa36702fdd7bfd43a60d7985ed811..fb75b03409990d8ad5c79e5e531de6b7917ce590 100644
--- a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecCherenkovAngles.h
+++ b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecCherenkovAngles.h
@@ -11,7 +11,7 @@ namespace Rich
     {
 
       /// Type for Cherenkov angles
-      using CherenkovAngles = Rich::Future::HypoData<double>;
+      using CherenkovAngles = Rich::Future::HypoData<float>;
 
       /// Cherenkov angles TES locations
       namespace CherenkovAnglesLocation
@@ -23,7 +23,7 @@ namespace Rich
       }
 
       /// type for expected Cherenkov resolutions
-      using CherenkovResolutions = Rich::Future::HypoData<double>;
+      using CherenkovResolutions = Rich::Future::HypoData<float>;
 
       /// Cherenkov resolutions TES locations
       namespace CherenkovResolutionsLocation
diff --git a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecGeomEfficiencies.h b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecGeomEfficiencies.h
index 70d1a85a3e7e3ab92314b3f9c79d51678857477e..5d784b3d4efc765777bdef64f1dff348f315d7e9 100644
--- a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecGeomEfficiencies.h
+++ b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecGeomEfficiencies.h
@@ -27,7 +27,7 @@ namespace Rich
       //=============================================================================
 
       /// Type for geometrical efficiencies
-      using GeomEffs = Rich::Future::HypoData<double>;
+      using GeomEffs = Rich::Future::HypoData<float>;
 
       /// geometrical efficiencies TES locations
       namespace GeomEffsLocation
@@ -39,7 +39,7 @@ namespace Rich
       //=============================================================================
 
       /// Type for the per PD Geom Effs for each mass hypothesis
-      using GeomEffsPerPD = Rich::ParticleArray< Rich::Map<LHCb::RichSmartID,double> >;
+      using GeomEffsPerPD = Rich::ParticleArray< Rich::Map<LHCb::RichSmartID,float> >;
 
       /// Vector of GeomEffsPerPD
       using GeomEffsPerPDVector = LHCb::STL::Vector<GeomEffsPerPD>;
diff --git a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecPhotonYields.h b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecPhotonYields.h
index f2212405108816a1d307f745c3ab23f9e315ab00..3a16f4430dd0caf8edcdd19a76776cf7125724f0 100644
--- a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecPhotonYields.h
+++ b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecPhotonYields.h
@@ -13,7 +13,7 @@ namespace Rich
     {
 
       /// Type for photon yield data
-      using PhotonYields = Rich::Future::HypoData<double>;
+      using PhotonYields = Rich::Future::HypoData<float>;
 
       /// photon yield TES locations
       namespace PhotonYieldsLocation
@@ -30,7 +30,7 @@ namespace Rich
       constexpr unsigned int NPhotonSpectraBins = 5;
 
       /// Type for photon spectra data
-      using PhotonSpectra = Rich::PhotonSpectra<double,NPhotonSpectraBins>;
+      using PhotonSpectra = Rich::PhotonSpectra<float,NPhotonSpectraBins>;
 
       /// photon spectra TES locations
       namespace PhotonSpectraLocation
diff --git a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecSIMDPixels.h b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecSIMDPixels.h
index be99a12b62e8d8cf83d64dd69fb162e76792aa65..aeea8b9962443e4d44980e8924de36c2d02b2d19 100644
--- a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecSIMDPixels.h
+++ b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecSIMDPixels.h
@@ -57,6 +57,8 @@ namespace Rich
         using SmartIDs = Rich::SIMD::STDArray<LHCb::RichSmartID,SIMDFP>;
         /// Type for index to original scalar cluster
         using ScIndex  = Rich::SIMD::Int32;
+        /// Selection mask
+        using Mask     = SIMDFP::MaskType;
 
       public:
 
@@ -67,14 +69,16 @@ namespace Rich
                    const Point &              gPos,
                    const Point &              lPos,
                    const SIMDFP &             effArea,
-                   const ScIndex &            scClusIn ) :
-          m_rich     ( rich     ),
-          m_side     ( side     ),
-          m_gloPos   ( gPos     ), 
-          m_locPos   ( lPos     ),
-          m_effArea  ( effArea  ),
-          m_smartID  ( smartIDs ),
-          m_scClusIn ( scClusIn ) { }
+                   const ScIndex &            scClusIn,
+                   const Mask &               mask ) :
+          m_rich      ( rich     ),
+          m_side      ( side     ),
+          m_gloPos    ( gPos     ), 
+          m_locPos    ( lPos     ),
+          m_effArea   ( effArea  ),
+          m_smartID   ( smartIDs ),
+          m_scClusIn  ( scClusIn ),
+          m_validMask ( mask     ) { }
 
       public:
 
@@ -99,6 +103,9 @@ namespace Rich
         /// Access the scalar cluster indices
         inline const ScIndex&  scClusIndex() const noexcept { return m_scClusIn; }
 
+        /// Access the validity mask
+        inline const Mask&     validMask() const noexcept { return m_validMask; }
+
       public:
 
         /// Implement ostream << method
@@ -138,6 +145,9 @@ namespace Rich
         /// Indices to the original scalar clusters
         ScIndex m_scClusIn = ScIndex(-1);
 
+        /// validity mask
+        Mask m_validMask { false };
+
       };
 
       /** @class SIMDPixelSummaries RichFutureRecEvent/RichRecSIMDPixels.h
diff --git a/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp b/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp
index 18f0c458d5debf509cc1d51e5261e674c0b43aae..fcaee5a16c40d929cbd54a94658a99587b9796ab 100644
--- a/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp
+++ b/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp
@@ -33,6 +33,7 @@ void SIMDPixelSummariesCreator::fill( const Rich::PDPixelCluster::Vector& cluste
   SIMDPixel::SIMDFP lX(0), lY(0), lZ(0); // local positions
   SIMDPixel::SIMDFP effArea(0);          // effective area
   SIMDPixel::ScIndex scClusIn(-1);       // indices to original clusters
+  SIMDPixel::Mask mask { false };        // selection mask
 
   // SIMD index
   std::size_t index = 0;
@@ -46,12 +47,27 @@ void SIMDPixelSummariesCreator::fill( const Rich::PDPixelCluster::Vector& cluste
       summaries.emplace_back( m_lastRich, m_lastSide, smartIDs,
                               SIMDPixel::Point(gX,gY,gZ),
                               SIMDPixel::Point(lX,lY,lZ),
-                              effArea, scClusIn ); 
+                              effArea, scClusIn, mask ); 
       // reset
       index = 0;
-      effArea = gX = gY = gZ = lX = lY = lZ = SIMDPixel::SIMDFP::Zero();
-      scClusIn = -1;
-      smartIDs.fill( LHCb::RichSmartID() );
+    };
+
+  // add padding info at the end of an incomplete SIMD pixel
+  auto addPadding = [&]()
+    {
+      // fill what has not been set with padding values
+      for ( std::size_t i = index; i < SIMDPixel::SIMDFP::Size; ++i )
+      {
+        effArea[i] = 0;
+        gX[i] = lX[i] = 0.0;
+        gY[i] = lY[i] = 0.0;
+        // Set default z to values roughly right for each RICH, to
+        // avoid precision issues in the photon reco later on.
+        gZ[i] = lZ[i] = ( Rich::Rich1 == m_lastRich ? 1600 : 10500 );
+        scClusIn[i] = -1;
+        mask[i] = false;
+        smartIDs[i] = LHCb::RichSmartID();
+      }
     };
 
   // Loop over all pixels
@@ -78,6 +94,8 @@ void SIMDPixelSummariesCreator::fill( const Rich::PDPixelCluster::Vector& cluste
     // If different RICH or side, save current data
     if ( UNLIKELY( rich != m_lastRich || side != m_lastSide ) )
     {
+      // add padding if needed
+      addPadding();
       // Save the current info
       savePix();
       // update RICH and panel
@@ -99,6 +117,8 @@ void SIMDPixelSummariesCreator::fill( const Rich::PDPixelCluster::Vector& cluste
     smartIDs[index] = cluster.primaryID();
     // scalar cluster index
     scClusIn[index] = clusIn;
+    // set mask OK for this pixel
+    mask[index] = true;
     
     // If this is the last index, push to container and start again
     if ( UNLIKELY( SIMDPixel::SIMDFP::Size-1 == index ) )
@@ -114,8 +134,14 @@ void SIMDPixelSummariesCreator::fill( const Rich::PDPixelCluster::Vector& cluste
  
   }
 
-  // handle padding values
-  if ( 0 != index ) { savePix(); }
+  // Save last one if needed
+  if ( 0 != index ) 
+  {
+    // add padding if needed
+    addPadding();
+    // Save the current info
+    savePix(); 
+  }
 
   // update ranges
   
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.cpp
index d528660076e5041bafeb5d6ed0cbfb1d5ace91bb..ae682ee90f325356ecf3f071f7583ce10108d55a 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.cpp
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.cpp
@@ -48,6 +48,10 @@ StatusCode BasePhotonReco::initialize()
     }
   }
 
+  _ri_debug << "Bias Corrs : " << m_ckBiasCorrs << endmsg;
+  
+  _ri_debug << "JO   Corrs : " << m_ckJOCorrs.value() << endmsg;
+
   // return
   return sc;
 }
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.h
index 89819395b0e743e0f8fd4384784b866799f8d233..26825e77d36cbf13a74b025c6af017724eb07235 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.h
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.h
@@ -191,27 +191,27 @@ namespace Rich
       protected: // Pre-sel options
   
         /// Min hit radius of interest around track centres for preselection
-        Gaudi::Property< RadiatorArray<double> > m_minROIPreSel
+        Gaudi::Property< RadiatorArray<float> > m_minROIPreSel
         { this, "PreSelMinTrackROI",  { 230,     0,       0 } };
         
         /// Max hit radius of interest around track centres for preselection
-        Gaudi::Property< RadiatorArray<double> > m_maxROIPreSel
+        Gaudi::Property< RadiatorArray<float> > m_maxROIPreSel
         { this,  "PreSelMaxTrackROI", { 540,    86,     165 } };
         
         /// N sigma for acceptance bands for preselection
-        Gaudi::Property< RadiatorArray<double> > m_nSigmaPreSel
+        Gaudi::Property< RadiatorArray<float> > m_nSigmaPreSel
         { this, "PreSelNSigma",       {  17,     6,      10 } };
 
         // cache values for speed
 
         /// Square of m_maxROI
-        RadiatorArray<double> m_maxROI2PreSel = {{}};
+        RadiatorArray<float> m_maxROI2PreSel = {{}};
         
         /// Square of m_minROI
-        RadiatorArray<double> m_minROI2PreSel = {{}};
+        RadiatorArray<float> m_minROI2PreSel = {{}};
         
         /// Internal cached parameter for speed
-        RadiatorArray<double> m_scalePreSel   = {{}};
+        RadiatorArray<float> m_scalePreSel   = {{}};
 
       protected: // reco options
 
@@ -220,39 +220,39 @@ namespace Rich
         { this, "CheckSideCrossing", { false, false, false } };
 
         /// Absolute minimum allowed Cherenkov Angle
-        Gaudi::Property< RadiatorArray<double> > m_minCKtheta
+        Gaudi::Property< RadiatorArray<float> > m_minCKtheta
         { this, "MinAllowedCherenkovTheta", { 0.150, 0.005, 0.005 },
             "The minimum allowed CK theta values for each radiator (Aero/R1Gas/R2Gas)" };
 
         /// Absolute maximum allowed Cherenkov Angle
-        Gaudi::Property< RadiatorArray<double> > m_maxCKtheta
+        Gaudi::Property< RadiatorArray<float> > m_maxCKtheta
         { this, "MaxAllowedCherenkovTheta", { 0.310, 0.075, 0.035 },
             "The maximum allowed CK theta values for each radiator (Aero/R1Gas/R2Gas)" };
         
         /// N sigma for acceptance bands
-        Gaudi::Property< RadiatorArray<double> > m_nSigma
+        Gaudi::Property< RadiatorArray<float> > m_nSigma
         { this, "NSigma", { 9.0,   5.25,  5.0   },
             "The CK theta # sigma selection range for each radiator (Aero/R1Gas/R2Gas)" };
 
         /** Cherenkov theta bias corrections, specific for each photon
          *  reconstruction method. */
-        RadiatorArray<double> m_ckBiasCorrs = {{}};
+        RadiatorArray<float> m_ckBiasCorrs = {{}};
 
       private:
 
         /** Job-Option Corrections applied to the reconstructed theta vales.
          *  By default 0. */
-        Gaudi::Property< RadiatorArray<double> > m_ckJOCorrs
+        Gaudi::Property< RadiatorArray<float> > m_ckJOCorrs
         { this, "CKThetaQuartzRefractCorrections", { 0.0, 0.0, 0.0 } };
 
         // Parameters for the scale factor calculations
 
         // The CK theta value
-        Gaudi::Property< RadiatorArray<double> > m_ckThetaScale
+        Gaudi::Property< RadiatorArray<float> > m_ckThetaScale
         { this, "ScaleFactorCKTheta", { 0.24, 0.052, 0.03 } };
 
         // The seperation the scale factors apply to
-        Gaudi::Property< RadiatorArray<double> > m_sepGScale
+        Gaudi::Property< RadiatorArray<float> > m_sepGScale
         { this, "ScaleFactorSepG", { 342,  75, 130 } };
 
       };
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichCKEstiFromRadiusPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichCKEstiFromRadiusPhotonReco.cpp
index 633ca0c6481bf9521cce16f39cb9c75d36089783..a1962ada9e743d4ed2c7c2714d22c33f006758dd 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichCKEstiFromRadiusPhotonReco.cpp
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichCKEstiFromRadiusPhotonReco.cpp
@@ -134,7 +134,7 @@ CKEstiFromRadiusPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& se
         const auto diff_y = ( segPSide.y() - pixPRad.y() );
 
         // compute the seperation squared
-        const auto sep2 = ( std::pow( diff_x, 2 ) + std::pow( diff_y, 2 ) );
+        const auto sep2 = ( ( diff_x * diff_x ) + ( diff_y * diff_y ) );
 
         // presel. default to rejected
         bool SelOK = false;
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.cpp
index b180b9759b79f53081f1a91344fe6c29582efc9a..5d88935806e5d4edf55b2862fb4a9bf19d99a4bf 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.cpp
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.cpp
@@ -12,12 +12,12 @@ CommonQuarticPhotonReco::CommonQuarticPhotonReco( const std::string& name,
 {
   // init
   m_rich.fill(nullptr);
-  m_ckBiasCorrs.fill(0);
   m_deBeam.fill(nullptr);
   
   // Corrections for the intrinsic biases
   //                  Aerogel      Rich1Gas   Rich2Gas
-  m_ckBiasCorrs = { -0.000358914, -7.505e-5, -4.287e-5 };
+  //m_ckBiasCorrs = { -0.000358914, -7.505e-5, -4.287e-5 }; // Run(1+2) HPD values
+  m_ckBiasCorrs = { 0.0,            2.0e-5,      -1.0e-5 }; // Run3 PMT values
 }
 
 //=============================================================================
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h
index 734ad5625b3cacedbefd90df1e29573cd5ca7018..80fc56cc92d3e6b319a9d2725ad586a4e484e8fa 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h
@@ -161,15 +161,15 @@ namespace Rich
         { this, "CheckPrimaryMirrorSegments", { false, false, false } };
 
         /// Minimum active segment fraction in each radiator
-        Gaudi::Property< RadiatorArray<double> > m_minActiveFrac
-        { this, "MinActiveFraction", { 0.2, 0.2, 0.2 } };
+        Gaudi::Property< RadiatorArray<float> > m_minActiveFrac
+        { this, "MinActiveFraction", { 0.2f, 0.2f, 0.2f } };
 
         /// Minimum tolerance for mirror reflection point during iterations
-        Gaudi::Property< RadiatorArray<double> > m_minSphMirrTolIt
+        Gaudi::Property< RadiatorArray<float> > m_minSphMirrTolIt
         { this, "MinSphMirrTolIt", 
-          { std::pow(0.10*Gaudi::Units::mm,2), 
-            std::pow(0.08*Gaudi::Units::mm,2), 
-            std::pow(0.05*Gaudi::Units::mm,2) } };
+          { std::pow( (float)(0.10f*Gaudi::Units::mm), 2 ), 
+            std::pow( (float)(0.08f*Gaudi::Units::mm), 2 ), 
+            std::pow( (float)(0.05f*Gaudi::Units::mm), 2 ) } };
 
         /// Flag to control if the secondary mirrors are treated as if they are completely flat
         Gaudi::Property< DetectorArray<bool> > m_treatSecMirrsFlat
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp
index 922576d1d7e24628ffd6e9e1d19f39e276cc7948..d5653ec65ca810e11dd5bda43990b1f0fd209a09 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp
@@ -65,7 +65,7 @@ QuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segments,
   
   // local position corrector
   // longer term need to remove this
-  const Rich::Rec::RadPositionCorrector<double> corrector;
+  const Rich::Rec::RadPositionCorrector<FP> corrector;
 
   // global photon index
   int photonIndex(-1);
@@ -145,8 +145,9 @@ QuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segments,
         const auto& segPanelPnt = tkLocPtn.point(side);
 
         // compute the seperation squared
-        const auto sep2 = ( std::pow( (pixP.x()-segPanelPnt.x()), 2 ) +
-                            std::pow( (pixP.y()-segPanelPnt.y()), 2 ) );
+        const auto dx   = pixP.x() - segPanelPnt.x();
+        const auto dy   = pixP.y() - segPanelPnt.y();
+        const auto sep2 = ( (dx*dx) + (dy*dy) );
 
         // presel. default to rejected
         bool SelOK = false;
@@ -380,11 +381,11 @@ QuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segments,
             virtDetPoint = gloPos - 2.0 * distance * secSegment->centreNormal();
           
             // solve the quartic using the new data
-            m_quarticSolver.solve<double,3,4>( emissionPoint,
-                                               sphSegment->centreOfCurvature(),
-                                               virtDetPoint,
-                                               sphSegment->radius(),
-                                               sphReflPoint );
+            m_quarticSolver.solve<FP,3,4>( emissionPoint,
+                                           sphSegment->centreOfCurvature(),
+                                           virtDetPoint,
+                                           sphSegment->radius(),
+                                           sphReflPoint );
           
             // (re)find the spherical mirror segment
             // CRJ - Is this needed ?
@@ -439,11 +440,11 @@ QuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segments,
               virtDetPoint = gloPos - ( normV * ( 2.0 * distance / normV.Mag2() ) );
 
               // solve the quartic using the new data
-              m_quarticSolver.solve<double,3,4>( emissionPoint,
-                                                 sphSegment->centreOfCurvature(),
-                                                 virtDetPoint,
-                                                 sphSegment->radius(),
-                                                 sphReflPoint );
+              m_quarticSolver.solve<FP,3,4>( emissionPoint,
+                                             sphSegment->centreOfCurvature(),
+                                             virtDetPoint,
+                                             sphSegment->radius(),
+                                             sphReflPoint );
 
               // (re)find the spherical mirror segment
               sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint );
@@ -531,7 +532,7 @@ QuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segments,
         // VDT and -O3 optimisation level. 
         // --------------------------------------------------------------------------------------
         photonDirection = (sphReflPoint-emissionPoint).Unit();
-        double thetaCerenkov(0), phiCerenkov(0);
+        FP thetaCerenkov(0), phiCerenkov(0);
         segment.angleToDirection( photonDirection, thetaCerenkov, phiCerenkov );
         //_ri_verbo << std::setprecision(9) 
         //          << " -> Phot Dir " << photonDirection << endmsg;
@@ -607,7 +608,7 @@ getBestGasEmissionPoint( const Rich::RadiatorType radiator,
                          Gaudi::XYZPoint & emissionPoint,
                          float & fraction ) const
 {
-  double alongTkFrac = 0.5;
+  FP alongTkFrac = 0.5;
   
   if ( radiator == Rich::Rich1Gas )
   {
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.h
index 8d08a17ade3b9edfef825a2115ec3a73a603ba4b..b6019c1cdaef1e8a2039f3974fd4eabf2cc6fdbb 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.h
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.h
@@ -70,6 +70,11 @@ namespace Rich
                             const SpacePointVector& localPoints,
                             const Relations::TrackToSegments::Vector& tkToSegRels ) const override;
 
+      private:
+
+        /// Floating point type to work with
+        using FP = double;
+
       private:
 
         /// Compute the best emission point for the gas radiators using
@@ -94,11 +99,11 @@ namespace Rich
                           Gaudi::XYZPoint & secReflPoint ) const
         {
           // solve quartic equation with nominal values and find spherical mirror reflection point
-          m_quarticSolver.solve<double,3,3>( emissionPoint,
-                                             m_rich[rich]->nominalCentreOfCurvature(side),
-                                             virtDetPoint,
-                                             m_rich[rich]->sphMirrorRadius(),
-                                             sphReflPoint );
+          m_quarticSolver.solve<FP,3,3>( emissionPoint,
+                                         m_rich[rich]->nominalCentreOfCurvature(side),
+                                         virtDetPoint,
+                                         m_rich[rich]->sphMirrorRadius(),
+                                         sphReflPoint );
           // find the spherical mirror segment
           //_ri_verbo << "   -> sphReflPoint " << sphReflPoint << endmsg;
           sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint );
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp
index 1438f4b108f59d6224d99fe62e75b8635b56230e..dc67af5f2b3b4c88898fdcf6033b3b3527990f1a 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp
@@ -66,7 +66,7 @@ StatusCode SIMDPhotonPredictedPixelSignal::initialize()
   }
 
   // lookup table interpolator for exp function
-  m_expFunc.init( m_minArg, 0, []( const double x ) { return std::exp(x); } );
+  //m_expFunc.init( m_minArg, 0, []( const double x ) { return std::exp(x); } );
 
   // cache min exp argument
   m_minArgSIMD = SIMDFP( m_minArg );
@@ -133,26 +133,28 @@ SIMDPhotonPredictedPixelSignal::operator()( const SIMDPixelSummaries& pixels,
     {
       if ( id != Rich::BelowThreshold )
       {
-        if ( tkCkAngles[id] >= 0.000001 )
+        // get the expect CK theta angle for this hypo
+        const auto tkA = tkCkAngles[id];
+        //if ( tkA >= 0.000001 )
         {
 
-          // Expected Cherenkov theta angle resolution
-          const SIMDFP thetaExpRes( tkCkRes[id] );
+          // 1 / Expected Cherenkov theta angle resolution
+          const SIMDFP thetaExpResInv ( SIMDFP::One() / SIMDFP(tkCkRes[id]) );
 
           // Expected CK theta for this hypo
-          const SIMDFP thetaExp( tkCkAngles[id] );
+          const SIMDFP thetaExp( tkA );
 
           // Theta diff / resolution
-          const SIMDFP thetaDiffOvRes = ( thetaReco - thetaExp ) / thetaExpRes;
+          const SIMDFP thetaDiffOvRes = ( thetaReco - thetaExp ) * thetaExpResInv;
             
           // compute the signal probability for this hypo
           // See note LHCB/98-040 page 11 equation 18
 
           // First the argument to exp() function
-          SIMDFP arg = SIMDFP(-0.5) * thetaDiffOvRes * thetaDiffOvRes;
+          SIMDFP arg = SIMDFP(-0.5f) * thetaDiffOvRes * thetaDiffOvRes;
 
           // selection mask
-          auto mask = arg > m_minExpArgF[rad] && phot.validityMask();
+          auto mask = phot.validityMask() && ( arg > m_minExpArgF[rad] );
 
           // If any are OK continue
           if ( any_of(mask) )
@@ -164,15 +166,21 @@ SIMDPhotonPredictedPixelSignal::operator()( const SIMDPixelSummaries& pixels,
             // compute exp(arg)
             const SIMDFP expArg = myexp(arg);
 
+            // Expected yield
+            const SIMDFP tkY( tkYields[id] );
+
             // The signal to set
-            SIMDFP & sig = sigs[id];
+            //SIMDFP & sig = sigs[id];
 
-            //const SIMDFP expArg = myexp(arg);
-            sig = AInd * expArg * SIMDFP(tkYields[id]) / thetaExpRes;
+            // Compute the signal
+            SIMDFP sig = AInd * expArg * tkY * thetaExpResInv;
             
             // Check min prob value
             mask &= sig > m_minPhotonProbSIMD[rad];
             sig.setZeroInverted(mask);
+
+            // save
+            sigs[id] = sig;
             
           } // min arg
 
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.h
index 652356fc01446ce569e4e63ca3a9a86c706d55c0..c66e0766925732f4b5c4c8bccc18f86b45f998d3 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.h
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.h
@@ -26,7 +26,7 @@
 // Rich Utils
 #include "RichUtils/ZipRange.h"
 #include "RichUtils/FastMaths.h"
-#include "RichUtils/LookupTableInterpolator.h"
+//#include "RichUtils/LookupTableInterpolator.h"
 
 // Detector Description
 #include "RichDet/DeRich1.h"
@@ -96,10 +96,10 @@ namespace Rich
         {
           // Fast VDT like function
           return Rich::Maths::fast_exp(x);
-          // Look up table
-          //return m_expFunc( x );
-          // testing ...
-          //return SIMDFP::One();
+          // Approximation
+          //return Rich::Maths::Approx::approx_exp(x);
+          // Even faster approximation
+          //return Rich::Maths::Approx::vapprox_exp(x);
         }
         
       private:
@@ -134,7 +134,7 @@ namespace Rich
         { this, "ScaleFactor", { 4.0, 4.0 } };
 
         /// Interpolator for exp function        
-        LookupTableInterpolator<float,5000> m_expFunc;
+        //LookupTableInterpolator<float,5000> m_expFunc;
         
       };
 
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
index 9b0fe61f2944048abe0bda5ada76868487362f4c..df0931746e76b51be831c358f07c0c1b37df41c9 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
@@ -75,7 +75,7 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
   auto & photons   = std::get<0>(outData);
   auto & relations = std::get<1>(outData);
   // guess at reserve size
-  const auto resSize = segments.size() * pixels.size() / 20;
+  const auto resSize = segments.size() * pixels.size() / 8;
   photons.reserve( resSize );
   relations.reserve( resSize );
   
@@ -83,6 +83,9 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
   // longer term need to remove this
   const Rich::Rec::RadPositionCorrector<SIMDPixel::SIMDFP> corrector;
 
+  // best secondary and spherical mirror data caches
+  MirrorData sphMirrors, secMirrors;
+
   // global photon index
   int photonIndex(-1);
   
@@ -117,8 +120,7 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
          SIMDPixel::Point(tkLocPtn.point(Rich::bottom)) }}; // also R2 right
 
       //_ri_verbo << endmsg;
-      //_ri_verbo << "Segment P=" << segment.bestMomentum() 
-      //          << " " << rich << " " << rad << endmsg;
+      //_ri_verbo << "Segment P=" << segment.bestMomentum() << " " << rich << " " << rad << endmsg;
 
       // get the best pixel range for this segment, based on where hits are expected
       const auto pixR = ( segFlags.inBothPanels() ? pixels.range(rich) :
@@ -144,6 +146,7 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
         //_ri_verbo << " -> Pixel IDs    " << pix.smartID() << endmsg;
         //_ri_verbo << " ->       GloPos " << pix.gloPos() << endmsg;
         //_ri_verbo << "          LocPos " << pix.locPos() << endmsg;
+        //_ri_verbo << "          Mask   " << pix.validMask() << endmsg;
 
         // side for this pack of pixels
         const auto side = pix.side();
@@ -157,14 +160,17 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
         ///           Probably correct but should check. 
         const auto& segPanelPnt = segSidePtns[side];
 
+        // Start mask off as the pixel mask
+        auto pixmask = pix.validMask();
+        
         // compute the seperation squared
         auto dx = pixP.x() - segPanelPnt.x();
         auto dy = pixP.y() - segPanelPnt.y();
         const auto sep2 = ( (dx*dx) + (dy*dy) );
 
         // Check overall boundaries
-        auto pixmask = ( sep2 < m_maxROI2PreSelSIMD[rad] && 
-                         sep2 > m_minROI2PreSelSIMD[rad] );
+        pixmask &= ( sep2 < m_maxROI2PreSelSIMD[rad] && 
+                     sep2 > m_minROI2PreSelSIMD[rad] );
         //_ri_verbo << "  -> sep2 = " << sep2 << " OK=" << pixmask << endmsg;
         if ( any_of(pixmask) )
         {
@@ -207,10 +213,6 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
         // fraction of segment path length accessible to the photon
         SIMDFP fraction{ SIMDFP::One() };
 
-        // best sec and spherical mirror data
-        MirrorData sphMirrors{{}};
-        MirrorData secMirrors{{}};
-
         // flag to say if these photon candidates are un-ambiguous - default to false
         SIMDFP::mask_type unambigPhoton { false };
 
@@ -261,7 +263,7 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
           findMirrorData( rich, side, virtDetPoint, emissionPoint2,
                           sphMirrors2, secMirrors2, 
                           sphReflPoint2, secReflPoint2, pixmask );
-          _ri_verbo << "  -> unambig2 " << pixmask << endmsg;
+          //_ri_verbo << "  -> unambig2 " << pixmask << endmsg;
           if ( UNLIKELY( none_of(pixmask) ) )
           {
             //_ri_debug << rad << " : Failed to reconstruct photon for end of segment" << endmsg;
@@ -371,6 +373,13 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
         }
         // --------------------------------------------------------------------------------------
 
+        // --------------------------------------------------------------------------------------
+        // make sure the mirror caches are up to date at this point
+        // --------------------------------------------------------------------------------------
+        secMirrors.update();
+        sphMirrors.update();
+        // --------------------------------------------------------------------------------------
+
         // --------------------------------------------------------------------------------------
         // Finally reconstruct the photon using best emission point and the best mirror segments
         // --------------------------------------------------------------------------------------
@@ -399,6 +408,8 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
             
             // (re)find the secondary mirror
             secMirrors.mirrors = m_mirrorSegFinder.get()->findSecMirror( rich, side, secReflPoint );
+            // update the cache
+            secMirrors.update();
 
             // Compute the virtual reflection point
  
@@ -416,23 +427,33 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
                                     ( normV.Z() * ( pix.gloPos().Z() - secReflPoint.Z() ) ) );
             virtDetPoint = pix.gloPos() - ( normV * ( SIMDFP(2.0) * distance / normV.Mag2() ) );
 
-            // solve the quartic using the new data ( try 2,3 )
-            m_quarticSolver.solve<SIMDFP,3,4>( emissionPoint,
+            // solve the quartic using the new data ( was 3,4 - try 2,3 )
+            m_quarticSolver.solve<SIMDFP,2,3>( emissionPoint,
                                                sphMirrors.getCoCs(),
                                                virtDetPoint,
                                                sphMirrors.getRoCs(),
                                                sphReflPoint );
 
-            // (re)find the spherical mirror segments
-            sphMirrors.mirrors = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint );
-            
             // for iterations after the first, see if we are still moving
             // If not, abort iterations early
             if ( iIt > m_nMinQits[rad] &&
-                 all_of( (last_mirror_point-secReflPoint).Mag2() < m_minSphMirrTolItSIMD[rad] ) ) { break; }
-
-            // store last sec mirror point
-            last_mirror_point = secReflPoint;
+                 all_of( (last_mirror_point-secReflPoint).Mag2() < m_minSphMirrTolItSIMD[rad] ) ) 
+            {
+              //_ri_verbo << "   -> Not Moving -> Abort iterations."<< endmsg;
+              break;
+            }
+
+            // Are we going to iterate again ?
+            if ( iIt+1 < m_nMaxQits[rad] )
+            {
+              //  (re)find the spherical mirror segments
+              //_ri_verbo << "   -> Update for sphMirrors using " << sphReflPoint << endmsg;
+              sphMirrors.mirrors = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint );
+              // update the caches
+              sphMirrors.update();
+              // store last sec mirror point
+              last_mirror_point = secReflPoint;
+            }
 
             // increment iteration counter
             ++iIt;
@@ -509,10 +530,9 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
         photonDirection = (sphReflPoint-emissionPoint).Unit();
         SIMDFP thetaCerenkov(SIMDFP::Zero()), phiCerenkov(SIMDFP::Zero());
         segment.angleToDirection( photonDirection, thetaCerenkov, phiCerenkov );
-        //_ri_verbo << std::setprecision(9) 
-        //          << " -> Phot Dir " << photonDirection << endmsg;
-        //_ri_verbo << std::setprecision(9) 
-        //          << "    theta,phi " << thetaCerenkov << " " << phiCerenkov << endmsg;
+        //_ri_verbo << std::setprecision(9) << " -> Phot Dir " << photonDirection << endmsg;
+        //_ri_verbo << std::setprecision(9) << "    theta    " << thetaCerenkov << endmsg;
+        //_ri_verbo << std::setprecision(9) << "    phi      " << phiCerenkov << endmsg;
         // --------------------------------------------------------------------------------------
 
         // --------------------------------------------------------------------------------------
@@ -566,11 +586,55 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
       
     } // segment loop
   } // track loop
+
+  // if ( photons.size() > resSize )
+  // {
+  //   warning() << "Photons reserve size too small. Scale = " 
+  //             << (float)( segments.size() * pixels.size() ) / photons.size() 
+  //             << endmsg;
+  // }
   
   // return the final data
   return outData;
 }
 
+//=============================================================================
+// Find mirror segments and reflection points for given data
+//=============================================================================
+void SIMDQuarticPhotonReco::findMirrorData( const Rich::DetectorType rich,
+                                            const Rich::Side side,
+                                            const SIMD::Point<FP> & virtDetPoint,
+                                            const SIMD::Point<FP> & emissionPoint,
+                                            Mirrors& primMirr,
+                                            Mirrors& secMirr,
+                                            SIMD::Point<FP> & sphReflPoint,
+                                            SIMD::Point<FP> & secReflPoint,
+                                            SIMDFP::mask_type & OK ) const
+{
+  // solve quartic equation with nominal values and find spherical mirror reflection point
+  // ( was 3,3 - try 2,2 )
+  m_quarticSolver.solve<SIMDFP,2,2>( emissionPoint,
+                                     m_rich[rich]->nominalCentreOfCurvatureSIMD(side),
+                                     virtDetPoint,
+                                     m_rich[rich]->sphMirrorRadiusSIMD(),
+                                     sphReflPoint );
+  // find the spherical mirror segments
+  //_ri_verbo << "   -> sphReflPoint " << sphReflPoint << endmsg;
+  primMirr = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint );
+  // Search for the secondary segment
+  // Direction vector from primary mirror point to virtual detection point
+  const auto dir ( virtDetPoint - sphReflPoint );
+  // find the sec mirror intersction point and secondary mirror segment
+  OK &= intersectPlane( sphReflPoint, dir,
+                        m_rich[rich]->nominalPlaneSIMD(side),
+                        secReflPoint );
+  if ( any_of(OK) )
+  {
+    // find the secondary mirrors
+    secMirr = m_mirrorSegFinder.get()->findSecMirror( rich, side, secReflPoint );
+  }
+}
+
 //=============================================================================
 // Compute the best emission point for the gas radiators using
 // the given spherical mirror reflection points
@@ -586,7 +650,7 @@ SIMDQuarticPhotonReco::getBestGasEmissionPoint( const Rich::RadiatorType radiato
 {
 
   // Default point for emission point is the middle
-  SIMDFP alongTkFrac = 0.5;
+  SIMDFP alongTkFrac( 0.5 );
 
   // Default to all OK
   SIMDFP::mask_type ok{true};
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.h
index cb006deed2792472455be93285a9faddcca3a732..337e43117cc770f0daa856df8c4ae80121b0b674 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.h
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.h
@@ -75,9 +75,11 @@ namespace Rich
 
       private:
 
-        // short cut to SIMD FP types
+        /// Basic precision (float)
         using FP      = SIMDCherenkovPhoton::FP;
+        /// SIMD version of FP
         using SIMDFP  = SIMDCherenkovPhoton::SIMDFP;
+        /// Type for container of mirror pointers. Same size as SIMDFP.
         using Mirrors = SIMD::STDArray< const DeRichSphMirror *, SIMDFP >;
 
       private:
@@ -85,45 +87,71 @@ namespace Rich
         /// Mirror Data Store
         class MirrorData
         {
+        public:
+          MirrorData() = default;
         public:
           /// Array of mirror DetElem pointers
           Mirrors mirrors{{}};
+        private: // data caches
+          /// The mirrors for which the caches are valid
+          Mirrors cache_mirrors{{}};
+          /// RoC values
+          SIMDFP          m_RoCs{};
+          /// CoC points
+          SIMD::Point<FP> m_CoCs;
+          /// Normal Planes
+          SIMD::Plane<FP> m_NormPs;
         public:
-          /// Get the SIMD RoC values for the current mirrors
-          SIMDFP getRoCs() const noexcept 
-          {
-            SIMDFP RoCs;
-            for ( std::size_t i = 0; i < SIMDFP::Size; ++i )
-            { RoCs[i] = mirrors[i]->radius(); }
-            return RoCs;
-          }
-          /// Get the SIMD centreNormalPlane for the current mirrors
-          SIMD::Plane<FP> getNormalPlane() const noexcept 
+          /// Update the cached information
+          void update() 
           {
-            SIMDFP A,B,C,D;
-            for ( std::size_t i = 0; i < SIMDFP::Size; ++i )
+            // Have the mirrors changed ?
+            if ( UNLIKELY( mirrors != cache_mirrors ) )
             {
-              const auto & p = mirrors[i]->centreNormalPlane();
-              A[i] = p.A();
-              B[i] = p.B();
-              C[i] = p.C();
-              D[i] = p.D();
+              SIMDFP A,B,C,D; // Plane params
+              SIMDFP X,Y,Z;   // CoC params
+              for ( std::size_t i = 0; i < SIMDFP::Size; ++i )
+              {
+                // the mirror pointer
+                const auto * m = mirrors[i];
+                // Update the RoC
+                m_RoCs[i] = m->radius();
+                // Update the plane params
+                const auto & p = m->centreNormalPlane();
+                A[i] = p.A();
+                B[i] = p.B();
+                C[i] = p.C();
+                D[i] = p.D();
+                // update the CoC params
+                const auto & n = m->centreOfCurvature();
+                X[i] = n.X();
+                Y[i] = n.Y();
+                Z[i] = n.Z();
+              }
+              // update the cache
+              cache_mirrors = mirrors;
+              m_NormPs      = SIMD::Plane<FP>(A,B,C,D);
+              m_CoCs        = SIMD::Point<FP>(X,Y,Z);
             }
-            return SIMD::Plane<FP>(A,B,C,D);
           }
-          /// Get the SIM CoC for the current mirrors
-          SIMD::Point<FP> getCoCs() const noexcept 
+        public:
+          ///< send to std::ostream
+          inline friend std::ostream& operator << ( std::ostream& s,
+                                                    const MirrorData& data )
           {
-            SIMDFP X,Y,Z;
-            for ( std::size_t i = 0; i < SIMDFP::Size; ++i )
-            {
-              const auto & p = mirrors[i]->centreOfCurvature();
-              X[i] = p.X();
-              Y[i] = p.Y();
-              Z[i] = p.Z();
-            }
-            return SIMD::Point<FP>(X,Y,Z);
+            return s << "{ Mirrors " << data.mirrors
+                     << " RoCs "  << data.getRoCs()
+                     << " CoCs "  << data.getCoCs()
+                     << " Plane " << data.getNormalPlane()
+                     << " }";
           }
+        public:
+          /// Get the SIMD RoC values for the current mirrors
+          const SIMDFP&                 getRoCs() const noexcept { return m_RoCs; }
+          /// Get the SIMD centreNormalPlane for the current mirrors
+          const SIMD::Plane<FP>& getNormalPlane() const noexcept { return m_NormPs; }
+          /// Get the SIM CoC for the current mirrors
+          const SIMD::Point<FP>&        getCoCs() const noexcept { return m_CoCs; }
         };
 
       private:
@@ -139,40 +167,15 @@ namespace Rich
                                                    SIMDFP & fraction ) const;
 
         /// Find mirror segments and reflection points for given data
-        inline void
-          findMirrorData( const Rich::DetectorType rich,
-                          const Rich::Side side,
-                          const SIMD::Point<FP> & virtDetPoint,
-                          const SIMD::Point<FP> & emissionPoint,
-                          Mirrors& primMirr,
-                          Mirrors& secMirr,
-                          SIMD::Point<FP> & sphReflPoint,
-                          SIMD::Point<FP> & secReflPoint,
-                          SIMDFP::mask_type & OK ) const
-        {
-          // solve quartic equation with nominal values and find spherical mirror reflection point
-          // ( try 2,2 )
-          m_quarticSolver.solve<SIMDFP,3,3>( emissionPoint,
-                                             m_rich[rich]->nominalCentreOfCurvatureSIMD(side),
-                                             virtDetPoint,
-                                             m_rich[rich]->sphMirrorRadiusSIMD(),
-                                             sphReflPoint );
-          // find the spherical mirror segments
-          //_ri_verbo << "   -> sphReflPoint " << sphReflPoint << endmsg;
-          primMirr = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint );
-          // Search for the secondary segment
-          // Direction vector from primary mirror point to virtual detection point
-          const auto dir ( virtDetPoint - sphReflPoint );
-          // find the sec mirror intersction point and secondary mirror segment
-          OK &= intersectPlane( sphReflPoint, dir,
-                                m_rich[rich]->nominalPlaneSIMD(side),
-                                secReflPoint );
-          if ( any_of(OK) )
-          {
-            // find the secondary mirrors
-            secMirr = m_mirrorSegFinder.get()->findSecMirror( rich, side, secReflPoint );
-          }
-        }
+        void findMirrorData( const Rich::DetectorType rich,
+                             const Rich::Side side,
+                             const SIMD::Point<FP> & virtDetPoint,
+                             const SIMD::Point<FP> & emissionPoint,
+                             Mirrors& primMirr,
+                             Mirrors& secMirr,
+                             SIMD::Point<FP> & sphReflPoint,
+                             SIMD::Point<FP> & secReflPoint,
+                             SIMDFP::mask_type & OK ) const;
 
         /// Check the final Cherenkov angles
         template < typename HYPODATA, typename FTYPE,
diff --git a/Rich/RichFutureRecPixelAlgorithms/src/RichSIMDPixelBackgroundsEstiAvHPD.cpp b/Rich/RichFutureRecPixelAlgorithms/src/RichSIMDPixelBackgroundsEstiAvHPD.cpp
index 7a8feb13b1a3fa43238a5fcdcfa251ce9bce8129..d89bb50ac1f07a5f1510d7390626153cce1cd88f 100644
--- a/Rich/RichFutureRecPixelAlgorithms/src/RichSIMDPixelBackgroundsEstiAvHPD.cpp
+++ b/Rich/RichFutureRecPixelAlgorithms/src/RichSIMDPixelBackgroundsEstiAvHPD.cpp
@@ -19,7 +19,7 @@ SIMDPixelBackgroundsEstiAvHPD( const std::string& name, ISvcLocator* pSvcLocator
                   KeyValue{ "GeomEffsPerPDLocation",         GeomEffsPerPDLocation::Default },
                   KeyValue{ "DetectablePhotonYieldLocation", PhotonYieldsLocation::Detectable },
                   KeyValue{ "RichSIMDPixelSummariesLocation", SIMDPixelSummariesLocation::Default } },
-                { KeyValue{ "PixelBackgroundsLocation",      PixelBackgroundsLocation::Default } } )
+                { KeyValue{ "PixelBackgroundsLocation",       SIMDPixelBackgroundsLocation::Default } } )
 {
   // debug
   //setProperty( "OutputLevel", MSG::VERBOSE );
@@ -78,6 +78,9 @@ SIMDPixelBackgroundsEstiAvHPD::operator()( const Relations::TrackToSegments::Vec
                                            const PhotonYields::Vector& detYieldsV,
                                            const SIMDPixelSummaries& pixels ) const
 {
+  // scalar type
+  using ScType = Rich::SIMD::DefaultScalarFP;
+
   // the backgrounds to return. Initialize to 0
   SIMDPixelBackgrounds          backgrounds( pixels.size(), SIMDFP::Zero()           );
 
@@ -201,13 +204,13 @@ SIMDPixelBackgroundsEstiAvHPD::operator()( const Relations::TrackToSegments::Vec
     
     int iter = 1;
     bool cont = true;
-    double rnorm = 0.0;
+    ScType rnorm = 0.0;
     while ( cont )
     {
       //_ri_debug << " -> Iteration " << iter << endmsg;
       
       int nBelow(0), nAbove(0);
-      double tBelow = 0.0;
+      ScType tBelow = 0.0;
       // loop over panels
       for ( auto& panelData : pdData[rich] )
       {
@@ -226,7 +229,7 @@ SIMDPixelBackgroundsEstiAvHPD::operator()( const Relations::TrackToSegments::Vec
               // First iteration, just set background for this HPD to the difference
               // between the observed and and expected number of hits in the HPD
               //_ri_debug << "  -> HPD " << pd << " obs. = " << obs << " exp. = " << exp << endmsg;
-              bkg = (double)iPD.obsSignal - iPD.expSignal;
+              bkg = (ScType)iPD.obsSignal - iPD.expSignal;
             }
             else
             {
@@ -258,7 +261,7 @@ SIMDPixelBackgroundsEstiAvHPD::operator()( const Relations::TrackToSegments::Vec
       {
         // we have some HPDs above and below expectation
         // calculate the amount of signal below per above HPD
-        rnorm = tBelow / ( static_cast<double>(nAbove) );
+        rnorm = tBelow / ( static_cast<ScType>(nAbove) );
         //_ri_debug << "   -> Correction factor per HPD above = " << rnorm << endmsg;
       }
       else
diff --git a/Rich/RichFutureRecSys/examples/IntelProfiler.py b/Rich/RichFutureRecSys/examples/IntelProfiler.py
new file mode 100644
index 0000000000000000000000000000000000000000..d87ec1752d8bd5fa8d4863b6ad6e69b9a1740f1a
--- /dev/null
+++ b/Rich/RichFutureRecSys/examples/IntelProfiler.py
@@ -0,0 +1,28 @@
+
+# Intel VTune auditor
+from Configurables import IntelProfilerAuditor
+profiler = IntelProfilerAuditor()
+#profiler.OutputLevel = DEBUG
+# We can skip some events
+profiler.StartFromEventN = 50
+#profiler.StopAtEventN = 500
+AuditorSvc().Auditors += [ profiler ]
+ApplicationMgr().AuditAlgorithms = True
+# Sequence which we need to profile. If empty, we profile all algorithms
+profiler.IncludeAlgorithms  = [ ]
+# Photon reco
+profiler.IncludeAlgorithms += [ "RichPhotonsDown",   "RichPhotonsUp",   "RichPhotonsLong" ]
+# ray tracing
+profiler.IncludeAlgorithms += [ "RichMassConesDown", "RichMassConesUp", "RichMassConesLong" ]
+# PID
+profiler.IncludeAlgorithms += [ "RichPIDDown",       "RichPIDUp",       "RichPIDLong" ]
+
+# Setup
+# source /cvmfs/projects.cern.ch/intelsw/psxe/linux/all-setup.sh
+
+# To avoid filling up /tmp
+# export TMPDIR=/group/rich/jonesc/IntelProfile/tmp
+
+# Example command line for the Intel profiler
+# amplxe-cl --collect general-exploration -start-paused -follow-child -target-duration-type=short -no-allow-multiple-runs -no-analyze-system -data-limit=1000 -q -- gaudirun.py -T  ~/LHCbCMake/Feature/Rec/Rich/RichFutureRecSys/examples/{RichFuture.py,data/MCXDstUpgradeFiles.py}
+
diff --git a/Rich/RichFutureRecSys/examples/RichFuture.py b/Rich/RichFutureRecSys/examples/RichFuture.py
index 8edc1f94b62458439d6763e0c1eee20779194c25..e1c09378588b6a9d6756a82785caa4d86d1cf6cc 100644
--- a/Rich/RichFutureRecSys/examples/RichFuture.py
+++ b/Rich/RichFutureRecSys/examples/RichFuture.py
@@ -1,7 +1,7 @@
 
 ###############################################################
 # Job options file
-#==============================================================
+###############################################################
 
 from Gaudi.Configuration import *
 from GaudiConfig.ControlFlow import seq
@@ -50,19 +50,19 @@ preloadGeom = True
 workAroundTrackTools = True
 
 # Input tracks
-#tkLocs = { "Long" : tkFilt.OutLongTracksLocation,
-#           "Down" : tkFilt.OutDownTracksLocation,
-#           "Up"   : tkFilt.OutUpTracksLocation }
+tkLocs = { "Long" : tkFilt.OutLongTracksLocation,
+           "Down" : tkFilt.OutDownTracksLocation,
+           "Up"   : tkFilt.OutUpTracksLocation }
 #tkLocs = { "All"  : "Rec/Track/Best" }
-tkLocs = { "Long" : tkFilt.OutLongTracksLocation }
+#tkLocs = { "Long" : tkFilt.OutLongTracksLocation }
 #tkLocs = { "Up" : tkFilt.OutUpTracksLocation }
 
 # Output PID
-#pidLocs =  { "Long" : "Rec/Rich/LongPIDs",
-#             "Down" : "Rec/Rich/DownPIDs",
-#             "Up"   : "Rec/Rich/UpPIDs" }
+pidLocs =  { "Long" : "Rec/Rich/LongPIDs",
+             "Down" : "Rec/Rich/DownPIDs",
+             "Up"   : "Rec/Rich/UpPIDs" }
 #pidLocs = { "All" : "Rec/Rich/PIDs" } 
-pidLocs = { "Long" : "Rec/Rich/LongPIDs" }
+#pidLocs = { "Long" : "Rec/Rich/LongPIDs" }
 #pidLocs = { "Up" : "Rec/Rich/UpPIDs" }
 
 # Merged output
@@ -124,13 +124,13 @@ all = GaudiSequencer( "All", Members = [ fetcher, odinDecode, richDecode,
 all.MeasureTime = True
 
 ApplicationMgr( TopAlg = [all],
-                EvtMax = 1000,      # events to be processed (default is 10)
+                EvtMax = 50000,      # events to be processed (default is 10)
                 #EvtSel = 'NONE', # do not use any event input
                 ExtSvc = ['ToolSvc', 'AuditorSvc' ],
                 AuditAlgorithms = True )
 
 EventSelector().PrintFreq = 100
-#LHCbApp().SkipEvents = 5
+#LHCbApp().SkipEvents = 5416
 
 CondDB() # Just to initialise
 LHCbApp()
@@ -153,6 +153,7 @@ RootHistCnv__PersSvc('RootHistCnv').ForceAlphaIds = True
 HistogramPersistencySvc().OutputFile = "RichFuture.root"
 
 # Auditors
-#AuditorSvc().Auditors += [ "FPEAuditor" ]
-#from Configurables import FPEAuditor
-#FPEAuditor().ActivateAt = ["Execute"]
+AuditorSvc().Auditors += [ "FPEAuditor" ]
+from Configurables import FPEAuditor
+FPEAuditor().ActivateAt = ["Execute"]
+#AuditorSvc().Auditors += [ "NameAuditor" ]
diff --git a/Rich/RichFutureRecSys/examples/RichFutureHive.py b/Rich/RichFutureRecSys/examples/RichFutureHive.py
index 8cf3311326f73fc9ab2b33cb26c1b53aa56669f7..e42a39701ae4bfbab3fca92ea1bb716861cec043 100644
--- a/Rich/RichFutureRecSys/examples/RichFutureHive.py
+++ b/Rich/RichFutureRecSys/examples/RichFutureHive.py
@@ -12,15 +12,15 @@ from GaudiConfig.ControlFlow import seq
 
 LHCbApp(EnableHive=True)
 
-from Configurables import AlgResourcePool, HiveWhiteBoard, ForwardSchedulerSvc
+#from Configurables import AlgResourcePool, HiveWhiteBoard, ForwardSchedulerSvc
 #AlgResourcePool().OutputLevel = VERBOSE
 
-scheduler = ForwardSchedulerSvc()
+#scheduler = ForwardSchedulerSvc()
 #scheduler.OutputLevel = VERBOSE
-whiteboard = HiveWhiteBoard("EventDataSvc")
+#whiteboard = HiveWhiteBoard("EventDataSvc")
 
-scheduler.MaxAlgosInFlight = 20
-scheduler.ThreadPoolSize   = 30
-whiteboard.EventSlots      = 6
+#scheduler.MaxAlgosInFlight = 20
+#scheduler.ThreadPoolSize   = 30
+#whiteboard.EventSlots      = 6
 
 #==============================================================
diff --git a/Rich/RichFutureRecSys/examples/data/MCXDstUpgradeFiles.py b/Rich/RichFutureRecSys/examples/data/MCXDstUpgradeFiles.py
index 1f0f435938c38e8173e0122a6e651ec900998e10..693325f656ba6edfb32fb10b71cffccab42cfd7b 100644
--- a/Rich/RichFutureRecSys/examples/data/MCXDstUpgradeFiles.py
+++ b/Rich/RichFutureRecSys/examples/data/MCXDstUpgradeFiles.py
@@ -4,8 +4,9 @@ from GaudiConf import IOHelper
 
 # Check what is available
 searchPaths = [
-    "/usera/jonesc/NFS/data/MC/Upgrade/XDST/",  # Cambridge
-    "/home/chris/LHCb/Data/MC/Upgrade/XDST/"    # CRJ's CernVM
+    "/usera/jonesc/NFS/data/MC/Upgrade/XDST/",         # Cambridge
+    "/home/chris/LHCb/Data/MC/Upgrade/XDST/",          # CRJ's CernVM
+    "/group/rich/jonesc/LHCb/Data/MC/Upgrade/XDST/"    # the pit
     ]
 
 print "Data Files :-"
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichDetailedTrSegMakerFromTracks.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichDetailedTrSegMakerFromTracks.cpp
index dd4fbffacae4d8b47a9ac6a4537eb6456e608317..37129f603a73a124c8f1a4e67b1749d550329fef 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichDetailedTrSegMakerFromTracks.cpp
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichDetailedTrSegMakerFromTracks.cpp
@@ -84,6 +84,11 @@ StatusCode DetailedTrSegMakerFromTracks::initialize()
   m_radiators.push_back( usedRads(Rich::Rich2Gas) ?
                          getDet<DeRichRadiator>(DeRichLocations::Rich2Gas) : nullptr );
 
+  for ( const auto rad : Rich::radiators() )
+  {
+    m_minRadLengthSq[rad] = std::pow( m_minRadLength[rad], 2 );
+  }
+
   if ( m_useStateProvider )
   {
     _ri_debug << "Will use StateProvider instead of extrapolator to move states" << endmsg;
@@ -495,7 +500,7 @@ constructSegments( const LHCb::Track * track,
     //---------------------------------------------------------------------------------------------
     // Radiator path length cut
     //---------------------------------------------------------------------------------------------
-    if ( (exitPState->position()-entryPState->position()).Mag2() < std::pow(m_minRadLength[rad],2) )
+    if ( (exitPState->position()-entryPState->position()).Mag2() < m_minRadLengthSq[rad] )
     {
       _ri_verbo << "    --> Path length too short -> rejecting segment" << endmsg;
       continue;
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichDetailedTrSegMakerFromTracks.h b/Rich/RichFutureRecTrackAlgorithms/src/RichDetailedTrSegMakerFromTracks.h
index 33e4e61d3c62105f1b653321e93d7664265a5709..f238bb99f7bade0dc663789e04ecab8d3350a0d2 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichDetailedTrSegMakerFromTracks.h
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichDetailedTrSegMakerFromTracks.h
@@ -216,6 +216,9 @@ namespace Rich
         { this, "MinRadiatorPathLength",
           { 0*Gaudi::Units::mm, 500*Gaudi::Units::mm, 1500*Gaudi::Units::mm } };
 
+        /// Cache min radius^2 at exit
+        RadiatorArray<double> m_minRadLengthSq = {{}};
+
         /// Flag to turn on/off the use of the TrackStateProvider to create missing states
         Gaudi::Property<bool> m_createMissingStates { this, "CreateMissingStates", true };
 
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichDetectablePhotonYields.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichDetectablePhotonYields.cpp
index 656172b58997055267fde6b54f414fe9c639858c..19065e5902b83b8ed9fa6d852aaae923cba56353 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichDetectablePhotonYields.cpp
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichDetectablePhotonYields.cpp
@@ -70,6 +70,9 @@ DetectablePhotonYields::operator()( const LHCb::RichTrackSegment::Vector& segmen
                                     const PhotonSpectra::Vector& emittedSpectra,
                                     const MassHypoRingsVector& massRings ) const
 {
+  // Scalar type to work with
+  using ScType = PhotonYields::Type;
+
   // make the data to return
   OutData data;
   auto & yieldV   = std::get<PhotonYields::Vector>(data);
@@ -110,7 +113,7 @@ DetectablePhotonYields::operator()( const LHCb::RichTrackSegment::Vector& segmen
     for ( const auto id : activeParticles() )
     {
       // the signal
-      double signal = 0;
+      ScType signal = 0;
 
       // Skip rings with no hits
       //_ri_debug << "  -> " << id << " #Rings=" << rings[id].size() << endmsg;
@@ -172,16 +175,16 @@ DetectablePhotonYields::operator()( const LHCb::RichTrackSegment::Vector& segmen
             const auto energy = emitSpectra.binEnergy(iEnBin);
 
             // Get weighted average PD Q.E.
-            double pdQEEff(0);
+            ScType pdQEEff(0);
             if ( !pdCount.empty() )
             {
               for ( const auto& PD : pdCount )
               {
                 // add up Q.E. eff
-                pdQEEff += (double)(PD.second) * (*((PD.first)->pdQuantumEff()))[energy*Gaudi::Units::eV]/100 ;
+                pdQEEff += (ScType)(PD.second) * (*((PD.first)->pdQuantumEff()))[energy*Gaudi::Units::eV]/100 ;
               }
               // normalise the result
-              pdQEEff /= (double)(totalInPD);
+              pdQEEff /= (ScType)(totalInPD);
             }
             else
             {
@@ -191,17 +194,17 @@ DetectablePhotonYields::operator()( const LHCb::RichTrackSegment::Vector& segmen
             //_ri_debug << std::setprecision(9) << "   -> pdQEEff : En=" << energy << " Eff=" << pdQEEff << endmsg;
 
             // Weighted primary mirror reflectivity
-            double primMirrRefl(0);
+            ScType primMirrRefl(0);
             if ( !primMirrCount.empty() )
             {
               for ( const auto& PM : primMirrCount )
               {
                 // add up mirror refl.
                 primMirrRefl +=
-                  (double)(PM.second) * (*(PM.first->reflectivity()))[energy*Gaudi::Units::eV];
+                  (ScType)(PM.second) * (*(PM.first->reflectivity()))[energy*Gaudi::Units::eV];
               }
               // normalise the result
-              primMirrRefl /= (double)(totalInPD);
+              primMirrRefl /= (ScType)(totalInPD);
             }
             else
             {
@@ -212,17 +215,17 @@ DetectablePhotonYields::operator()( const LHCb::RichTrackSegment::Vector& segmen
             //_ri_debug << std::setprecision(9) << "   -> primMirrRefl : En=" << energy << " Eff=" << primMirrRefl << endmsg;
 
             // Weighted secondary mirror reflectivity
-            double secMirrRefl(0);
+            ScType secMirrRefl(0);
             if ( !secMirrCount.empty() )
             {
               for ( const auto& SM : secMirrCount )
               {
                 // add up mirror refl.
                 secMirrRefl +=
-                  (double)(SM.second) * (*(SM.first->reflectivity()))[energy*Gaudi::Units::eV];
+                  (ScType)(SM.second) * (*(SM.first->reflectivity()))[energy*Gaudi::Units::eV];
               }
               // normalise the result
-              secMirrRefl /= (double)(totalInPD);
+              secMirrRefl /= (ScType)(totalInPD);
             }
             else
             {
@@ -240,8 +243,8 @@ DetectablePhotonYields::operator()( const LHCb::RichTrackSegment::Vector& segmen
             
             // The Quartz window efficiency
             sig *= ( energy <= 0 ? 0 :
-                     vdt::fast_exp( -m_qWinZSize[rich] /
-                                    (*(m_riches[rich]->gasWinAbsLength()))[energy*Gaudi::Units::eV] ) );
+                     Rich::Maths::fast_exp( -m_qWinZSize[rich] /
+                                            (*(m_riches[rich]->gasWinAbsLength()))[energy*Gaudi::Units::eV] ) );
 
             //_ri_debug << std::setprecision(9) << "    -> Final Eff = " << totEff << " | DetSignal = " << sig << endmsg;
 
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichDetectablePhotonYields.h b/Rich/RichFutureRecTrackAlgorithms/src/RichDetectablePhotonYields.h
index 840672bc1bf396d6074c57cdbe5b9e38490bab78..9873b15b21ac49100bfe478a6b78acd317f5a0cd 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichDetectablePhotonYields.h
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichDetectablePhotonYields.h
@@ -24,6 +24,7 @@
 #include "RichUtils/ZipRange.h"
 #include "RichUtils/StlArray.h"
 #include "RichUtils/RichMap.h"
+#include "RichUtils/FastMaths.h"
 
 // RichDet
 #include "RichDet/DeRich1.h"
@@ -31,9 +32,6 @@
 #include "RichDet/DeRichPD.h"
 #include "RichDet/DeRichSphMirror.h"
 
-// VDT
-#include "vdt/exp.h"
-
 namespace Rich
 {
   namespace Future
@@ -94,10 +92,10 @@ namespace Rich
         DetectorArray<const DeRich*> m_riches = {{}};
 
         /// Cached value storing product of quartz window eff. and digitisation pedestal loss
-        double m_qEffPedLoss{0};
+        float m_qEffPedLoss{0};
 
         /// Thickness of windows in each RICH
-        DetectorArray<double> m_qWinZSize = {{}};
+        DetectorArray<float> m_qWinZSize = {{}};
 
       };
 
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichEmittedPhotonYields.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichEmittedPhotonYields.cpp
index 149d8b802b49e20254efa19959f09af8a26774e8..4471c2ed9cee652292db8961b91fddd92a85f26b 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichEmittedPhotonYields.cpp
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichEmittedPhotonYields.cpp
@@ -85,7 +85,7 @@ EmittedPhotonYields::operator()( const LHCb::RichTrackSegment::Vector& segments
     for ( const auto id : activeParticles() )
     {
       // the signal value
-      double signal = 0;
+      PhotonYields::Type signal = 0;
 
       // Skip below threshold
       if ( Rich::BelowThreshold != id )
@@ -114,7 +114,7 @@ EmittedPhotonYields::operator()( const LHCb::RichTrackSegment::Vector& segments
             const auto & radInts = segment.radIntersections();
             
             // normalise over each intersection
-            double totPlength(0), waveIndepTrans(0);
+            PhotonYields::Type totPlength(0), waveIndepTrans(0);
             if ( !radInts.empty() )
             {
               // average energy for this range
@@ -126,7 +126,7 @@ EmittedPhotonYields::operator()( const LHCb::RichTrackSegment::Vector& segments
                 const auto pLen = R.pathLength();
                 const auto absL = R.radiator()->absorption()->value(avEn);
                 totPlength       += pLen;
-                waveIndepTrans   += pLen * vdt::fast_exp( -pLen / absL );
+                waveIndepTrans   += pLen * Rich::Maths::fast_exp( -pLen / absL );
               }
               if ( totPlength>0 ) { waveIndepTrans /= totPlength; }
               
@@ -136,7 +136,7 @@ EmittedPhotonYields::operator()( const LHCb::RichTrackSegment::Vector& segments
           }
 
           // bin signal
-          const auto binSignal = length * ( nPhot > 0 ? nPhot : 0 );
+          const auto binSignal = length * ( nPhot > 0 ? nPhot : 0.0 );
           //_ri_verbo << std::setprecision(9) 
           //          << " bin " << iEnBin << " nPhot " << nPhot 
           //          << " binSignal " << binSignal << endmsg;
@@ -230,7 +230,7 @@ StatusCode EmittedPhotonYields::umsUpdate()
     const auto RH = RE02/RF;
     const auto RM = RE + (2.0*RC);
     const auto RS = RG + (2.0*RC);
-    const auto RT = sqrt( 0.25*RM*RM - RH*RS );
+    const auto RT = std::sqrt( 0.25*RM*RM - RH*RS );
     const auto RXSP = std::sqrt( (RM/2.0 + RT)/RH );
     const auto RXSM = std::sqrt( (RM/2.0 - RT)/RH );
     REP[rad] = std::sqrt(RE02) * RXSP;
@@ -248,9 +248,9 @@ StatusCode EmittedPhotonYields::umsUpdate()
       ( const Rich::RadiatorType rad, const double energy )
       {
         const auto A = ( RXSPscale[rad] *
-                         vdt::fast_log( (REP[rad]+energy)/(REP[rad]-energy) ) );
+                         Rich::Maths::fast_log( (REP[rad]+energy)/(REP[rad]-energy) ) );
         const auto B = ( RXSMscale[rad] *
-                         vdt::fast_log( (REM[rad]+energy)/(REM[rad]-energy) ) );
+                         Rich::Maths::fast_log( (REM[rad]+energy)/(REM[rad]-energy) ) );
         return X[rad] * ( A - B );
       } );
 
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichEmittedPhotonYields.h b/Rich/RichFutureRecTrackAlgorithms/src/RichEmittedPhotonYields.h
index 073efb3d4542c5f3c21c5b98f88a9bcc58503c32..8b767a8fbb4c42101013d341d05454c09ecd5800 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichEmittedPhotonYields.h
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichEmittedPhotonYields.h
@@ -16,6 +16,7 @@
 
 // Utils
 #include "RichUtils/RichTrackSegment.h"
+#include "RichUtils/FastMaths.h"
 
 // Interfaces
 #include "RichInterfaces/IRichDetParameters.h"
@@ -24,10 +25,6 @@
 #include "RichDet/DeRich1.h"
 #include "RichDet/DeRichAerogelRadiator.h"
 
-// VDT
-#include "vdt/exp.h"
-#include "vdt/log.h"
-
 namespace Rich
 {
   namespace Future
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.cpp
index 2dde84e00df38cc0bb15b58579d7234a65b6ccf9..2282d4ffc53b56a16d0e01c182d8e3c144cfbcc3 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.cpp
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.cpp
@@ -25,6 +25,9 @@ GeomEffCKMassRings::operator()( const LHCb::RichTrackSegment::Vector& segments,
                                 const CherenkovAngles::Vector& ckAngles,
                                 const MassHypoRingsVector& massRings ) const
 {
+  // Scalar type
+  using ScType = GeomEffs::Type;
+
   // make the data to return
   OutData data;
   auto & geomEffsV      = std::get<GeomEffs::Vector>(data);
@@ -60,19 +63,19 @@ GeomEffCKMassRings::operator()( const LHCb::RichTrackSegment::Vector& segments,
     for ( const auto id : activeParticles() )
     {
       //The efficiency
-      double eff = 0;
+      ScType eff = 0;
       
       // above threshold ?
       if ( ckAngles[id] > 0 )
       {
         // count the number of detected points
-        int nDetect(0);
+        std::size_t nDetect(0);
 
         // number of points
         const auto nPoints = massRings[id].size();
 
         // PD increment
-        const auto pdInc = 1 / static_cast<float>(nPoints);
+        const auto pdInc = 1.0 / static_cast<float>(nPoints);
 
         // Loop over the points of the mass ring
         for ( const auto& P : massRings[id] )
@@ -95,7 +98,7 @@ GeomEffCKMassRings::operator()( const LHCb::RichTrackSegment::Vector& segments,
         }
 
         // compute the final eff
-        eff = static_cast<float>(nDetect)/static_cast<float>(nPoints);
+        eff = static_cast<ScType>(nDetect)/static_cast<ScType>(nPoints);
 
       }
 
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp
index dcacb7aaf86a6f96fc63e897758b6e2b2887411e..ee12a248d449efd035daabd065a2a76ed863dcac 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp
@@ -99,6 +99,9 @@ RayTraceCherenkovCones::operator()( const LHCb::RichTrackSegment::Vector& segmen
   // longer term need to remove this
   const Rich::Rec::RadPositionCorrector<double> corrector;
 
+  // Comparison result for good ray tracings
+  const RayTracingUtils::SIMDResult::Results goodRes( LHCb::RichTraceMode::InPDPanel );
+
   // loop over the input data
   for ( const auto && data : Ranges::ConstZip(segments,ckAngles) )
   {
@@ -153,22 +156,26 @@ RayTraceCherenkovCones::operator()( const LHCb::RichTrackSegment::Vector& segmen
         unsigned int nPhot(0), nOK(0);
         for ( const auto && data : Ranges::ConstZip(results,m_cosSinPhiV[rad]) )
         {
-          const auto & res    = std::get<0>(data);
+          // results
+          const auto & res = std::get<0>(data);
+
+          // Count total photons
+          nPhot += SIMDFP::Size;
+          // Count good photons
+          nOK   += ( res.result >= goodRes ).count();
+          // bailout check
+          if ( 0 == nOK && nPhot >= m_nBailout[rad] ) { break; }  
+
+          // angles
           const auto & cosphi = std::get<1>(data);
 
           // detection point (SIMD)
           const auto & gPSIMD = res.detectionPoint;
 
-          // break out flag
-          bool goOn{true};
-
           // Scalar loop
           for ( std::size_t i = 0; i < SIMDFP::Size; ++i )
           {
             
-            // count photons
-            ++nPhot;
-
             // Scalar hit point
             const Gaudi::XYZPoint gP( gPSIMD.x()[i], gPSIMD.y()[i], gPSIMD.z()[i] ); 
 
@@ -185,18 +192,9 @@ RayTraceCherenkovCones::operator()( const LHCb::RichTrackSegment::Vector& segmen
             //_ri_verbo << std::setprecision(4) 
             //          << LHCb::RichTraceMode::RayTraceResult(res.result[i]) << " " << gP << endmsg;
             //_ri_verbo << rings[id].back() << endmsg;
-            
-            // count raytraces that are in HPD panel
-            if ( res.result[i] >= LHCb::RichTraceMode::InPDPanel ) { ++nOK; }
-            
-            // bailout check
-            if ( 0 == nOK && nPhot >= m_nBailout[rad] ) { goOn = false; break; }  
 
           } // scalar loop
 
-          // break if need be
-          if ( UNLIKELY(!goOn) ) { break; }
-
         } // results loop
         
         // if no good hits empty the container
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackCherenkovAngles.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackCherenkovAngles.cpp
index 39392f7bba00bb13bf6367f83431f725b6e0048f..6768db665f9101533316d9c0992abb11bc1e52b3 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackCherenkovAngles.cpp
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackCherenkovAngles.cpp
@@ -28,6 +28,7 @@ TrackCherenkovAnglesBase::run( const LHCb::RichTrackSegment::Vector& segments,
                                const PhotonSpectra::Vector& tkSpectra,
                                const PhotonYields::Vector& tkYields ) const
 {
+
   // make the data to return
   CherenkovAngles::Vector anglesV;
   anglesV.reserve( segments.size() );
@@ -50,7 +51,7 @@ TrackCherenkovAnglesBase::run( const LHCb::RichTrackSegment::Vector& segments,
     {
 
       // the angle
-      double angle = 0;
+      CherenkovAngles::Type angle = 0;
 
       // Is this hypo above threshold ?
       if ( yields[id] > 0 )
@@ -73,7 +74,7 @@ TrackCherenkovAnglesBase::run( const LHCb::RichTrackSegment::Vector& segments,
             //          << " bin " << iEnBin << " " << beta << " " << refIn << endmsg;
             if ( temp>1 ) 
             {
-              const auto f  = vdt::fast_acos(1.0/temp);
+              const auto f  = Rich::Maths::fast_acos(1.0/temp);
               const auto en = (spectra.energyDist(id))[iEnBin];
               //_ri_verbo << std::setprecision(9) << "     " << " " << f << " " << en << endmsg;
               angle += ( en * f );
@@ -85,7 +86,7 @@ TrackCherenkovAnglesBase::run( const LHCb::RichTrackSegment::Vector& segments,
         }
 
       }
-      
+ 
       // save the final angle for this hypo
       //_ri_verbo << std::setprecision(9) << id << " CK theta " << angle << endmsg;
       angles.setData(id,angle);
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackCherenkovAngles.h b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackCherenkovAngles.h
index 41bdb6981bc86dcd822daf255ae0fe785a629860..74172b99c54ca32becc37a0cfdb4d3a9ae8ec125 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackCherenkovAngles.h
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackCherenkovAngles.h
@@ -6,9 +6,6 @@
 #include <algorithm>
 #include <iomanip>
 
-// VDT math
-#include "vdt/asin.h"
-
 // Gaudi Functional
 #include "GaudiAlg/Transformer.h"
 
@@ -22,6 +19,7 @@
 // Utils
 #include "RichUtils/RichTrackSegment.h"
 #include "RichUtils/ZipRange.h"
+#include "RichUtils/FastMaths.h"
 
 // Interfaces
 #include "RichInterfaces/IRichRefractiveIndex.h"
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.cpp
index fbca806dcfcf0b93f1384cd6ef2066f2e5162982..33df2d18b0f0b8450b861981d89287bae842c58d 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.cpp
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.cpp
@@ -102,6 +102,8 @@ operator()( const LHCb::RichTrackSegment::Vector& segments,
 
     // momentum for this segment, in GeV units
     const auto ptot = segment.bestMomentumMag() / GeV;
+    // 1 / momentum
+    const auto ptotInv = 1.0 / ptot;
 
     // check if segment is valid
     if ( aboveThreshold && ptot > 0 )
@@ -138,8 +140,9 @@ operator()( const LHCb::RichTrackSegment::Vector& segments,
       if ( effL > 0 )
       {
         const auto multScattCoeff = ( m_scatt * std::sqrt(effL) *
-                                      ( 1.0 + 0.038*vdt::fast_log(effL) ) );
-        const auto scattErr = 2.0 * std::pow(multScattCoeff/ptot,2);
+                                      ( 1.0 + 0.038 * Rich::Maths::fast_log(effL) ) );
+        const auto scatCOvP = multScattCoeff * ptotInv;
+        const auto scattErr = 2.0 * scatCOvP * scatCOvP;
         //_ri_debug << std::setprecision(9) << " Scattering error = " << scattErr << endmsg;
         res2 += scattErr;
       }
@@ -175,7 +178,7 @@ operator()( const LHCb::RichTrackSegment::Vector& segments,
         // Uses the supplied reference contribution for the given pixel area, and scales
         // according to the area for the associated PD  
         //-------------------------------------------------------------------------------
-        auto pdErr = std::pow( m_pdErr[rad], 2 );
+        auto pdErr = m_pdErr[rad] * m_pdErr[rad];
         // Are we using the full treatment to deal with PD size difference ?
         if ( UNLIKELY(m_fullPDAreaTreatment[rad]) )
         {
@@ -215,7 +218,9 @@ operator()( const LHCb::RichTrackSegment::Vector& segments,
         {
           
           // cache tan(cktheta)
-          const auto tanCkExp = vdt::fast_tan(angles[hypo]);
+          const auto tanCkExp = Rich::Maths::fast_tan(angles[hypo]);
+          // 1 /  tan(cktheta)
+          const auto tanCkExpInv = 1.0 / tanCkExp;
           //_ri_debug << std::setprecision(9) << "  " << hypo
           //          << " ckExp = " << angles[hypo]
           //          << " tanCkExp = " << tanCkExp << endmsg;
@@ -223,12 +228,13 @@ operator()( const LHCb::RichTrackSegment::Vector& segments,
           //-------------------------------------------------------------------------------
           // chromatic error
           //-------------------------------------------------------------------------------
-          const auto index      = m_refIndex.get()->refractiveIndex(segment.radIntersections(),
-                                                                    segment.avPhotonEnergy());
-          const auto chromFact  = ( index>0 ? m_detParams.get()->refIndexSD(rad)/index : 0 );
-          const auto chromatErr = std::pow( chromFact / tanCkExp, 2 );
+          const auto index        = m_refIndex.get()->refractiveIndex(segment.radIntersections(),
+                                                                      segment.avPhotonEnergy());
+          const auto chromFact    = ( index>0 ? m_detParams.get()->refIndexSD(rad)/index : 0 );
+          const auto chromatErr   = chromFact * tanCkExpInv;
+          const auto chromatErrSq = chromatErr * chromatErr;
           //_ri_debug << std::setprecision(9) << "         Chromatic err = " << chromatErr << endmsg;
-          hypo_res2 += chromatErr;
+          hypo_res2 += chromatErrSq;
           //-------------------------------------------------------------------------------
 
           //-------------------------------------------------------------------------------
@@ -236,18 +242,18 @@ operator()( const LHCb::RichTrackSegment::Vector& segments,
           //-------------------------------------------------------------------------------
           const auto mass2      = richPartProps()->massSq(hypo)/(GeV*GeV);
           const auto massFactor = mass2 / ( mass2 + (ptot*ptot) );
-          const auto momErr     = ( segment.entryErrors().errP2()/(GeV*GeV) *
-                                    std::pow( massFactor / ( ptot * tanCkExp ), 2 ) );
-          //_ri_debug << std::setprecision(9) << "         momentum err = " << momErr << endmsg;
-          hypo_res2 += momErr;
+          const auto A          = massFactor * ptotInv * tanCkExpInv;
+          const auto momErrSq   = ( (segment.entryErrors().errP2()/(GeV*GeV)) * A * A );
+          //_ri_debug << std::setprecision(9) << "         momentum err = " << momErrSq << endmsg;
+          hypo_res2 += momErrSq;
           //-------------------------------------------------------------------------------
 
           //-------------------------------------------------------------------------------
           // Global Scale factor
           //-------------------------------------------------------------------------------
-          hypo_res2 *= std::pow( m_scale[rad], 2 );
+          hypo_res2 *= m_scale[rad] * m_scale[rad];
           //-------------------------------------------------------------------------------
-          
+
         } // theta > 0
 
         // Compute the final resolution
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.h b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.h
index 65c12c19442dfe6f1bae219ee55b5786ca7a642f..e086350b4f1ee52609c19d04c0c95ca73b4132e9 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.h
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.h
@@ -26,6 +26,7 @@
 #include "RichUtils/RichTrackSegment.h"
 #include "RichUtils/ZipRange.h"
 #include "RichUtils/RichGeomFunctions.h"
+#include "RichUtils/FastMaths.h"
 
 // Interfaces
 #include "RichInterfaces/IRichRefractiveIndex.h"
@@ -35,10 +36,6 @@
 #include "DetDesc/ITransportSvc.h"
 #include "DetDesc/TransportSvcException.h"
 
-// VDT
-#include "vdt/log.h"
-#include "vdt/tan.h"
-
 namespace Rich
 {
   namespace Future
diff --git a/Rich/RichRecTests/src/PhotonReco/main.icpp b/Rich/RichRecTests/src/PhotonReco/main.icpp
index a13821b494123339fb8b3c3fe76241c28c7fa5b3..670956f1848c8ad112695508e4c3fbf522630b6c 100644
--- a/Rich/RichRecTests/src/PhotonReco/main.icpp
+++ b/Rich/RichRecTests/src/PhotonReco/main.icpp
@@ -93,7 +93,7 @@ solveV( const typename DATA::Vector & dataV, POINT & sphReflPoint )
 
 int main ( int /*argc*/, char** /*argv*/ )
 {
-  const std::size_t nPhotons = 96e2; // must be devisible by 16
+  const std::size_t nPhotons = 48e2; // must be devisible by 16
   
   Data<double>::Vector dataVD;
   Data<float>::Vector  dataVF;
diff --git a/Rich/RichRecTests/src/RayTracing/main.icpp b/Rich/RichRecTests/src/RayTracing/main.icpp
index 8356899652828de3cef3ba6ddcc118cc93ea4c4a..0d3729b9b31e70b37f88628efacc557b51152d6e 100644
--- a/Rich/RichRecTests/src/RayTracing/main.icpp
+++ b/Rich/RichRecTests/src/RayTracing/main.icpp
@@ -79,7 +79,7 @@ unsigned long long int __attribute__((noinline)) rayTrace( const DATA & in_dataV
 
   unsigned long long int best_dur{ 99999999999999999 };
 
-  const unsigned int nTests = 10000;
+  const unsigned int nTests = 50000;
   //std::cout << "Running " << nTests << " tests on " << in_dataV.size() << " data objects" << std::endl;
   DATA dataV; // working data copy;
   for ( unsigned int i = 0; i < nTests; ++i )
@@ -98,7 +98,7 @@ unsigned long long int __attribute__((noinline)) rayTrace( const DATA & in_dataV
 int main ( int /*argc*/, char** /*argv*/ )
 {
 
-  const unsigned int nPhotons = 96e2;
+  const unsigned int nPhotons = 48e2;
   //const unsigned int nPhotons = 8;
   std::cout << "Creating " << nPhotons << " random photons ..." << std::endl;
 
@@ -137,9 +137,10 @@ int main ( int /*argc*/, char** /*argv*/ )
     //std::cout << "VectorClass double " << resVClass << " speedup " << (float)resROOT/(float)resVClass << std::endl;
 
     // Vc float
-    Data<VcXYZPointD,VcXYZVectorD,VcPlane3DD,Vc::double_v>::Vector VC_dataV( nPhotons / Vc::double_v::Size );
+    Data<VcXYZPointD,VcXYZVectorD,VcPlane3DD,Vc::double_v>::Vector VC_dataV( nPhotons );
     auto resVC = rayTrace( VC_dataV );
-    std::cout << "Vc double          " << resVC << " speedup " << (float)resROOT/(float)resVC << std::endl;
+    std::cout << "Vc double          " << resVC << " per scalar speedup " 
+              << Vc::double_v::Size * (float)resROOT/(float)resVC << std::endl;
 
     // make sure we are not optimized away
     asm volatile ("" : "+x"(resROOT));
@@ -152,7 +153,7 @@ int main ( int /*argc*/, char** /*argv*/ )
   //if ( false )
   {
 
-   // Data in 'ROOT' SMatrix float format
+    // Data in 'ROOT' SMatrix float format
     Data<LHCbXYZPointF,LHCbXYZVectorF,LHCbPlane3DF,LHCbXYZPointF::Scalar>::Vector root_dataV( nPhotons );
     auto resROOT = rayTrace( root_dataV );
     std::cout << "ROOT float         " << resROOT << std::endl;
@@ -168,9 +169,10 @@ int main ( int /*argc*/, char** /*argv*/ )
     //std::cout << "VectorClass float  " << resVClass << " speedup " << (float)resROOT/(float)resVClass << std::endl;
 
     // Vc float
-    Data<VcXYZPointF,VcXYZVectorF,VcPlane3DF,Vc::float_v>::Vector VC_dataV( nPhotons / Vc::float_v::Size );
+    Data<VcXYZPointF,VcXYZVectorF,VcPlane3DF,Vc::float_v>::Vector VC_dataV( nPhotons );
     auto resVC = rayTrace( VC_dataV );
-    std::cout << "Vc float           " << resVC << " speedup " << (float)resROOT/(float)resVC << std::endl;
+    std::cout << "Vc float           " << resVC << " per scalar speedup " 
+              << Vc::float_v::Size * (float)resROOT/(float)resVC << std::endl;
 
     // make sure we are not optimized away
     asm volatile ("" : "+x"(resROOT));
diff --git a/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h b/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h
index b6195d77a600012eadb066e19b0660ad1121a25e..f8ffdc275626d4dff78b7e6e652c05f13d350385 100755
--- a/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h
+++ b/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h
@@ -314,8 +314,6 @@ namespace Rich
       template < class TYPE >
       inline TYPE f4( const TYPE (&a)[5], const TYPE& x ) const
       {
-        //return (a0 * x*x*x*x + a1 * x*x*x + a2 * x*x + a3 * x + a4);
-        //A bit more FMA friendly
         return ( ( ( ( ( ( ( a[0] * x ) + a[1] ) * x ) + a[2] ) * x ) + a[3] ) * x ) + a[4];
       }
 
@@ -327,6 +325,7 @@ namespace Rich
       inline void  
       evalPolyHorner( const TYPE (&a)[ORDER+1], TYPE (&res)[DIFFGRADE+1], const TYPE & x ) const
       {
+
         for ( std::size_t i = 0; i <= DIFFGRADE; ++i )
         {
           res[i] = a[0];
@@ -334,11 +333,16 @@ namespace Rich
 
         for ( std::size_t j = 1; j <= ORDER; ++j )
         {
-          res[0] = res[0] * x + a[j];
-          const int l = (ORDER - j) > DIFFGRADE ? DIFFGRADE : ORDER - j;
-          for ( int i = 1; i <= l ; ++i )
+          //res[0] = res[0] * x + a[j];
+          res[0] *= x;
+          res[0] += a[j];
+          //const std::size_t l = ( ORDER-j > DIFFGRADE ? DIFFGRADE : ORDER-j );
+          const auto l = std::min( ORDER-j, DIFFGRADE );
+          for ( std::size_t i = 1; i <= l ; ++i )
           {
-            res[i] = res[i] * x + res[i-1];
+            res[i] *= x;
+            res[i] += res[i-1];
+            //res[i] = res[i] * x + res[i-1];
           }
         }
 
@@ -346,8 +350,9 @@ namespace Rich
         TYPE l = 1.0;
         for ( std::size_t i = 2; i <= DIFFGRADE; ++i )
         {
-          l *= i;
-          res[i] = res[i] * l;
+          l      *= i;
+          res[i] *= l;
+          //res[i] = res[i] * l;
         }
 
       }