From fdea2cc16751150f59eefc2e99c82bb3f7f49a99 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Thu, 30 Nov 2017 14:55:46 +0000
Subject: [PATCH 01/29] Switch from look up table interpolater, which has cache
 miss problems, to a fast approximate implementation of log(exp(x)-1)

---
 .../RichSIMDGlobalPIDLikelihoodMinimiser.cpp  | 10 +-
 .../RichSIMDGlobalPIDLikelihoodMinimiser.h    | 97 +++++++++++++++++--
 2 files changed, 96 insertions(+), 11 deletions(-)

diff --git a/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.cpp b/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.cpp
index a953a060196..8da2f4afb9f 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 03708a0a5f9..85b3d3319ae 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
@@ -202,11 +278,20 @@ namespace Rich
                     typename std::enable_if<!std::is_arithmetic<TYPE>::value>::type * = nullptr >
           inline TYPE fast_logExp( const TYPE & x ) const noexcept
           {
+
+            //info() << "test " << x << " = " << approx_logExp(x) << " " << full_logExp(x) << endmsg;
+
             // 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 +399,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
-- 
GitLab


From 0849421ac1e4cea330c75591ba803099df0086fe Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Fri, 1 Dec 2017 14:08:33 +0000
Subject: [PATCH 02/29] Switch by default to the fastest log(exp(x)-1) option

---
 .../src/RichSIMDGlobalPIDLikelihoodMinimiser.h             | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.h b/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.h
index 85b3d3319ae..868f7092183 100644
--- a/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.h
+++ b/Rich/RichFutureGlobalPID/src/RichSIMDGlobalPIDLikelihoodMinimiser.h
@@ -269,8 +269,8 @@ namespace Rich
             // 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 );
+            //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
@@ -278,9 +278,6 @@ namespace Rich
                     typename std::enable_if<!std::is_arithmetic<TYPE>::value>::type * = nullptr >
           inline TYPE fast_logExp( const TYPE & x ) const noexcept
           {
-
-            //info() << "test " << x << " = " << approx_logExp(x) << " " << full_logExp(x) << endmsg;
-
             // Use the fast VDT inspired methods
             //return full_logExp(x);
 
-- 
GitLab


From b3a85e95176ef70e5a6493ef52498d4370cf65d2 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Fri, 1 Dec 2017 14:08:56 +0000
Subject: [PATCH 03/29] update mathematica notebook

---
 .../mathematica/Log(Exp(x)-1).nb              | 212 ++++++++++--------
 1 file changed, 115 insertions(+), 97 deletions(-)

diff --git a/Rich/RichFutureGlobalPID/mathematica/Log(Exp(x)-1).nb b/Rich/RichFutureGlobalPID/mathematica/Log(Exp(x)-1).nb
index c356ae46e2a..6fcd54dce91 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  ]]
 }
 ]
-- 
GitLab


From 605ccb35407bf2059fcd6aeef0c5c4b56fa45eec Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Fri, 1 Dec 2017 14:10:09 +0000
Subject: [PATCH 04/29] Add ROOT script to visualize the various approximate
 log(exp(x)-1) options

---
 Rich/RichFutureGlobalPID/scripts/drawfuncs.C | 88 ++++++++++++++++++++
 1 file changed, 88 insertions(+)
 create mode 100644 Rich/RichFutureGlobalPID/scripts/drawfuncs.C

diff --git a/Rich/RichFutureGlobalPID/scripts/drawfuncs.C b/Rich/RichFutureGlobalPID/scripts/drawfuncs.C
new file mode 100644
index 00000000000..dbe89ab38a0
--- /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");
+
+}
-- 
GitLab


From 1d68970874a1143e8696e382f8cc6e88ee37e255 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Fri, 1 Dec 2017 20:20:23 +0000
Subject: [PATCH 05/29] Add fast exp options

---
 .../src/RichSIMDPhotonPredictedPixelSignal.cpp       |  2 +-
 .../src/RichSIMDPhotonPredictedPixelSignal.h         | 12 ++++++------
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp
index 5694626383a..007ebd6e6d7 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 );
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.h
index 652356fc014..c66e0766925 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;
         
       };
 
-- 
GitLab


From c1211a3989add4e546ebaece7e03a7e0a28800b3 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sat, 2 Dec 2017 00:05:37 +0100
Subject: [PATCH 06/29] Convert some data to float

---
 .../RichFutureRecEvent/RichRecCherenkovAngles.h               | 4 ++--
 .../RichFutureRecEvent/RichRecGeomEfficiencies.h              | 4 ++--
 .../RichFutureRecEvent/RichRecPhotonYields.h                  | 4 ++--
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecCherenkovAngles.h b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecCherenkovAngles.h
index 1447eba3a28..fb75b034099 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 70d1a85a3e7..5d784b3d4ef 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 f2212405108..3a16f4430dd 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
-- 
GitLab


From 2ee7968f73e0b195af981c5d8d64bd4f8ee09f45 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sat, 2 Dec 2017 00:06:58 +0100
Subject: [PATCH 07/29] Use aliases to define scalar types to use

---
 .../src/RichDetectablePhotonYields.cpp        | 27 ++++++++++---------
 .../src/RichDetectablePhotonYields.h          |  8 +++---
 .../src/RichEmittedPhotonYields.cpp           | 14 +++++-----
 .../src/RichEmittedPhotonYields.h             |  5 +---
 .../src/RichGeomEffCKMassRings.cpp            |  7 +++--
 .../src/RichTrackCherenkovAngles.cpp          |  7 ++---
 .../src/RichTrackCherenkovAngles.h            |  4 +--
 7 files changed, 36 insertions(+), 36 deletions(-)

diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichDetectablePhotonYields.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichDetectablePhotonYields.cpp
index 076c552071c..fa7a1b2b37c 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 840672bc1bf..9873b15b21a 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 229d0f1c9e7..257b99d0b0b 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 073efb3d454..8b767a8fbb4 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 f000df9dc67..5a32d93ecc8 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,7 +63,7 @@ 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 )
@@ -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/RichTrackCherenkovAngles.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackCherenkovAngles.cpp
index 3367353f7f4..a84d116522a 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 41bdb6981bc..74172b99c54 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"
-- 
GitLab


From 1c754dc5a58e9a7d37c054ce6e37ccc95f0fd12b Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sat, 2 Dec 2017 00:07:33 +0100
Subject: [PATCH 08/29] Use aliases to define scalar types to use

---
 .../src/RichSIMDPixelBackgroundsEstiAvHPD.cpp       | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/Rich/RichFutureRecPixelAlgorithms/src/RichSIMDPixelBackgroundsEstiAvHPD.cpp b/Rich/RichFutureRecPixelAlgorithms/src/RichSIMDPixelBackgroundsEstiAvHPD.cpp
index a7b2b1eeb4c..f3560666cc8 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
-- 
GitLab


From 6a410cc4fb23ff12baab5bdc045e80df48a74e38 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sat, 2 Dec 2017 00:08:17 +0100
Subject: [PATCH 09/29] update examaple options

---
 Rich/RichFutureRecSys/examples/RichFuture.py         |  2 +-
 Rich/RichFutureRecSys/examples/RichFutureHive.py     | 12 ++++++------
 .../examples/data/MCXDstUpgradeFiles.py              |  3 ++-
 3 files changed, 9 insertions(+), 8 deletions(-)

diff --git a/Rich/RichFutureRecSys/examples/RichFuture.py b/Rich/RichFutureRecSys/examples/RichFuture.py
index 8edc1f94b62..1f3c2ca979c 100644
--- a/Rich/RichFutureRecSys/examples/RichFuture.py
+++ b/Rich/RichFutureRecSys/examples/RichFuture.py
@@ -124,7 +124,7 @@ all = GaudiSequencer( "All", Members = [ fetcher, odinDecode, richDecode,
 all.MeasureTime = True
 
 ApplicationMgr( TopAlg = [all],
-                EvtMax = 1000,      # events to be processed (default is 10)
+                EvtMax = 5000,      # events to be processed (default is 10)
                 #EvtSel = 'NONE', # do not use any event input
                 ExtSvc = ['ToolSvc', 'AuditorSvc' ],
                 AuditAlgorithms = True )
diff --git a/Rich/RichFutureRecSys/examples/RichFutureHive.py b/Rich/RichFutureRecSys/examples/RichFutureHive.py
index 8cf3311326f..e42a39701ae 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 1f0f435938c..7dff41f7b2c 100644
--- a/Rich/RichFutureRecSys/examples/data/MCXDstUpgradeFiles.py
+++ b/Rich/RichFutureRecSys/examples/data/MCXDstUpgradeFiles.py
@@ -5,7 +5,8 @@ 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
+    "/home/chris/LHCb/Data/MC/Upgrade/XDST/",   # CRJ's CernVM
+    "/group/rich/LHCb/Data/MC/Upgrade/XDST/"    # the pit
     ]
 
 print "Data Files :-"
-- 
GitLab


From a9a729007ecfd99513e8afdcaf74e5f5aed3bc4e Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sat, 2 Dec 2017 00:18:37 +0100
Subject: [PATCH 10/29] update the pixel signal algorithm

---
 .../RichSIMDPhotonPredictedPixelSignal.cpp    | 28 ++++++++++++-------
 1 file changed, 18 insertions(+), 10 deletions(-)

diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp
index 007ebd6e6d7..757366f0912 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDPhotonPredictedPixelSignal.cpp
@@ -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
 
-- 
GitLab


From 40d33aec01d9b9e65dfa2ca131a9107aa399f2e0 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sat, 2 Dec 2017 00:24:02 +0100
Subject: [PATCH 11/29] change data location at the pit

---
 Rich/RichFutureRecSys/examples/data/MCXDstUpgradeFiles.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/Rich/RichFutureRecSys/examples/data/MCXDstUpgradeFiles.py b/Rich/RichFutureRecSys/examples/data/MCXDstUpgradeFiles.py
index 7dff41f7b2c..693325f656b 100644
--- a/Rich/RichFutureRecSys/examples/data/MCXDstUpgradeFiles.py
+++ b/Rich/RichFutureRecSys/examples/data/MCXDstUpgradeFiles.py
@@ -4,9 +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
-    "/group/rich/LHCb/Data/MC/Upgrade/XDST/"    # the pit
+    "/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 :-"
-- 
GitLab


From e80434ea6ae257f9ee520e21a3242685f00de0c0 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sat, 2 Dec 2017 13:32:25 +0000
Subject: [PATCH 12/29] Avoid some pow calls

---
 .../src/RichCKEstiFromRadiusPhotonReco.cpp    |  2 +-
 .../src/RichQuarticPhotonReco.cpp             |  5 +--
 .../src/RichDetailedTrSegMakerFromTracks.cpp  |  7 +++-
 .../src/RichDetailedTrSegMakerFromTracks.h    |  3 ++
 ...ichTrackFunctionalCherenkovResolutions.cpp | 36 +++++++++++--------
 .../RichTrackFunctionalCherenkovResolutions.h |  5 +--
 6 files changed, 35 insertions(+), 23 deletions(-)

diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichCKEstiFromRadiusPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichCKEstiFromRadiusPhotonReco.cpp
index 8f152b1f372..29d98bf4b0a 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/RichQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp
index 5bc489330ed..4d260f93bc0 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp
@@ -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;
diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichDetailedTrSegMakerFromTracks.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichDetailedTrSegMakerFromTracks.cpp
index 0fef6712ac6..f95546e2e2a 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 33e4e61d3c6..f238bb99f7b 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/RichTrackFunctionalCherenkovResolutions.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.cpp
index 7ea7c69bd35..2398927af0a 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 65c12c19442..e086350b4f1 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
-- 
GitLab


From cbc8c015f798aafe6e2dbd893d66d1056d3b305c Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sat, 2 Dec 2017 18:05:38 +0000
Subject: [PATCH 13/29] tune newton solver

---
 Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h b/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h
index b6195d77a60..19a659efd11 100755
--- a/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h
+++ b/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h
@@ -335,9 +335,12 @@ 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 )
+          const auto ordermj  = ORDER - j;
+          const std::size_t l = ( ordermj > DIFFGRADE ? DIFFGRADE : ordermj );
+          for ( std::size_t i = 1; i <= l ; ++i )
           {
+            //res[i] *= x;
+            //res[i] += res[i-1];
             res[i] = res[i] * x + res[i-1];
           }
         }
@@ -347,6 +350,7 @@ namespace Rich
         for ( std::size_t i = 2; i <= DIFFGRADE; ++i )
         {
           l *= i;
+          //res[i] *= l;
           res[i] = res[i] * l;
         }
 
-- 
GitLab


From 9b1f2e148fbf4a9a5ddf74e53b019bf49bc7141e Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sat, 2 Dec 2017 20:00:26 +0100
Subject: [PATCH 14/29] tune the horner method

---
 .../RichRecUtils/QuarticSolverNewton.h        | 23 ++++++++++---------
 1 file changed, 12 insertions(+), 11 deletions(-)

diff --git a/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h b/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h
index 19a659efd11..f8ffdc27562 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,14 +333,16 @@ namespace Rich
 
         for ( std::size_t j = 1; j <= ORDER; ++j )
         {
-          res[0] = res[0] * x + a[j];
-          const auto ordermj  = ORDER - j;
-          const std::size_t l = ( ordermj > DIFFGRADE ? DIFFGRADE : ordermj );
+          //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] *= x;
-            //res[i] += res[i-1];
-            res[i] = res[i] * x + res[i-1];
+            res[i] *= x;
+            res[i] += res[i-1];
+            //res[i] = res[i] * x + res[i-1];
           }
         }
 
@@ -349,9 +350,9 @@ namespace Rich
         TYPE l = 1.0;
         for ( std::size_t i = 2; i <= DIFFGRADE; ++i )
         {
-          l *= i;
-          //res[i] *= l;
-          res[i] = res[i] * l;
+          l      *= i;
+          res[i] *= l;
+          //res[i] = res[i] * l;
         }
 
       }
-- 
GitLab


From 73ec5171fc322b806696c0018241955ff3ada8c6 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sat, 2 Dec 2017 20:00:47 +0100
Subject: [PATCH 15/29] Add some options for the interl profiler

---
 Rich/RichFutureRecSys/examples/RichFuture.py | 47 ++++++++++++++++----
 1 file changed, 38 insertions(+), 9 deletions(-)

diff --git a/Rich/RichFutureRecSys/examples/RichFuture.py b/Rich/RichFutureRecSys/examples/RichFuture.py
index 1f3c2ca979c..df2743ff95a 100644
--- a/Rich/RichFutureRecSys/examples/RichFuture.py
+++ b/Rich/RichFutureRecSys/examples/RichFuture.py
@@ -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,7 +124,7 @@ all = GaudiSequencer( "All", Members = [ fetcher, odinDecode, richDecode,
 all.MeasureTime = True
 
 ApplicationMgr( TopAlg = [all],
-                EvtMax = 5000,      # 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 )
@@ -156,3 +156,32 @@ HistogramPersistencySvc().OutputFile = "RichFuture.root"
 #AuditorSvc().Auditors += [ "FPEAuditor" ]
 #from Configurables import FPEAuditor
 #FPEAuditor().ActivateAt = ["Execute"]
+
+# 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}
+
+
-- 
GitLab


From 5d07cec7e515c38c27690327ef596ca806a61f21 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sat, 2 Dec 2017 22:37:10 +0000
Subject: [PATCH 16/29] count good ray traced photons outside scalar loop

---
 .../src/RichRayTraceCherenkovCones.cpp        | 30 +++++++++----------
 1 file changed, 14 insertions(+), 16 deletions(-)

diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp
index ca11245f36d..8c61d2313e5 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
-- 
GitLab


From 952d72a9869ce9615f2f3064d38bd8ee735259b2 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sun, 3 Dec 2017 14:15:03 +0000
Subject: [PATCH 17/29] Improvements to the photon reconstruction algorithm.
 Change default properties to float, reduce the number of iterations in the
 newton solver and implement a more efficient cache of the mirror parameters.

---
 .../src/RichBasePhotonReco.cpp                |   4 +
 .../src/RichBasePhotonReco.h                  |  26 ++---
 .../src/RichCommonQuarticPhotonReco.cpp       |   4 +-
 .../src/RichCommonQuarticPhotonReco.h         |   4 +-
 .../src/RichQuarticPhotonReco.cpp             |  26 ++---
 .../src/RichQuarticPhotonReco.h               |  15 ++-
 .../src/RichSIMDQuarticPhotonReco.cpp         |  35 ++++--
 .../src/RichSIMDQuarticPhotonReco.h           | 102 ++++++++++--------
 8 files changed, 127 insertions(+), 89 deletions(-)

diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.cpp
index d528660076e..ae682ee90f3 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 89819395b0e..26825e77d36 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/RichCommonQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.cpp
index b180b9759b7..5d88935806e 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 734ad5625b3..49c5b167d3d 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h
@@ -161,11 +161,11 @@ namespace Rich
         { this, "CheckPrimaryMirrorSegments", { false, false, false } };
 
         /// Minimum active segment fraction in each radiator
-        Gaudi::Property< RadiatorArray<double> > m_minActiveFrac
+        Gaudi::Property< RadiatorArray<float> > m_minActiveFrac
         { this, "MinActiveFraction", { 0.2, 0.2, 0.2 } };
 
         /// 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), 
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp
index 4d260f93bc0..29d6308d5e6 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);
@@ -381,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 ?
@@ -440,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 );
@@ -532,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;
@@ -608,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 8d08a17ade3..b6019c1cdae 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/RichSIMDQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
index db8514d2c0d..b382d640826 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
@@ -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);
   
@@ -207,10 +210,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 };
 
@@ -371,6 +370,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 +405,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 +424,28 @@ 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;
+            // Are we going to iterate again ?
+            if ( iIt+1 < m_nMaxQits[rad] )
+            {
+              //  (re)find the spherical mirror segments
+              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;
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.h
index cb006deed27..3071906125c 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,60 @@ 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();
+              // update the cache
+              cache_mirrors = mirrors;
+              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();
+              }
+              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);
           }
+        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
-          SIMD::Point<FP> getCoCs() const noexcept 
-          {
-            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);
-          }
+          const SIMD::Point<FP>&        getCoCs() const noexcept { return m_CoCs; }
         };
 
       private:
@@ -139,20 +156,19 @@ 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
+        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,
+          // ( was 3,3 - try 2,2 )
+          m_quarticSolver.solve<SIMDFP,2,2>( emissionPoint,
                                              m_rich[rich]->nominalCentreOfCurvatureSIMD(side),
                                              virtDetPoint,
                                              m_rich[rich]->sphMirrorRadiusSIMD(),
-- 
GitLab


From 81905f607c25a44bbf249e80dd4d8a63c8971bff Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sun, 3 Dec 2017 17:24:18 +0000
Subject: [PATCH 18/29] Move options for Intel profiler to seperate options
 file

---
 .../examples/IntelProfiler.py                 | 28 ++++++++++
 Rich/RichFutureRecSys/examples/RichFuture.py  | 51 ++++---------------
 2 files changed, 39 insertions(+), 40 deletions(-)
 create mode 100644 Rich/RichFutureRecSys/examples/IntelProfiler.py

diff --git a/Rich/RichFutureRecSys/examples/IntelProfiler.py b/Rich/RichFutureRecSys/examples/IntelProfiler.py
new file mode 100644
index 00000000000..d87ec1752d8
--- /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 df2743ff95a..7cd50ef8184 100644
--- a/Rich/RichFutureRecSys/examples/RichFuture.py
+++ b/Rich/RichFutureRecSys/examples/RichFuture.py
@@ -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
@@ -118,13 +118,13 @@ all = GaudiSequencer( "All", Members = [ fetcher, odinDecode, richDecode,
                                          tkUnpack, tkFilt, RichRec ] )
 
 # Uncomment to enable monitoring
-#all.Members += [ RichMoni ]
-#all.Members += [ mcPUnpack, mcRiSumUnP, RichCheck ]
+all.Members += [ RichMoni ]
+all.Members += [ mcPUnpack, mcRiSumUnP, RichCheck ]
 
 all.MeasureTime = True
 
 ApplicationMgr( TopAlg = [all],
-                EvtMax = 50000,      # events to be processed (default is 10)
+                EvtMax = 5000,      # events to be processed (default is 10)
                 #EvtSel = 'NONE', # do not use any event input
                 ExtSvc = ['ToolSvc', 'AuditorSvc' ],
                 AuditAlgorithms = True )
@@ -156,32 +156,3 @@ HistogramPersistencySvc().OutputFile = "RichFuture.root"
 #AuditorSvc().Auditors += [ "FPEAuditor" ]
 #from Configurables import FPEAuditor
 #FPEAuditor().ActivateAt = ["Execute"]
-
-# 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}
-
-
-- 
GitLab


From bc64f80157b946b677d5877bcd905184a0260fde Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Tue, 5 Dec 2017 08:44:59 +0000
Subject: [PATCH 19/29] Add pixel selection mask to SIMD pixels to identify
 padding entries. Make default z for padding positions 5000mm to help avoid
 precision issues later on.

---
 .../RichFutureRecEvent/RichRecSIMDPixels.h    | 26 +++++++++++++------
 .../src/RichRecSIMDPixels.cpp                 | 25 +++++++++++++-----
 2 files changed, 37 insertions(+), 14 deletions(-)

diff --git a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecSIMDPixels.h b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecSIMDPixels.h
index be99a12b62e..aeea8b99624 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 18f0c458d5d..aaae787467c 100644
--- a/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp
+++ b/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp
@@ -28,11 +28,19 @@ void SIMDPixelSummariesCreator::fill( const Rich::PDPixelCluster::Vector& cluste
   // last side
   Rich::Side         m_lastSide { Rich::InvalidSide     };
 
+  // Default z position for invalid entries. 
+  // NB : Not zero to help precision problems downstream in photon reco.
+  const SIMDPixel::FP defZ = 5000.0;
+
+  // Default SmartID to fill
+  const LHCb::RichSmartID defID();
+
   // Working Vc values
-  SIMDPixel::SIMDFP gX(0), gY(0), gZ(0); // global positions
-  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::SIMDFP gX(0), gY(0), gZ(defZ); // global positions
+  SIMDPixel::SIMDFP lX(0), lY(0), lZ(defZ); // 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,11 +54,14 @@ 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();
+      effArea = SIMDPixel::SIMDFP::Zero();
+      gX = gY = lX = lY = SIMDPixel::SIMDFP::Zero();
+      gZ = lZ = SIMDPixel::SIMDFP(defZ);
       scClusIn = -1;
+      mask = SIMDPixel::Mask{false};
       smartIDs.fill( LHCb::RichSmartID() );
     };
 
@@ -99,6 +110,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 ) )
-- 
GitLab


From 88b77a63d72cfded0f0cfdaee8d707b27cc65880 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Tue, 5 Dec 2017 08:46:33 +0000
Subject: [PATCH 20/29] Start photon mask as pixel mask.

---
 .../src/RichSIMDQuarticPhotonReco.cpp         | 59 ++++++++++++++++---
 .../src/RichSIMDQuarticPhotonReco.h           | 53 +++++++----------
 2 files changed, 72 insertions(+), 40 deletions(-)

diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
index b382d640826..f7be20440c4 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
@@ -120,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) :
@@ -147,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();
@@ -160,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) )
         {
@@ -260,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;
@@ -430,16 +433,21 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
                                                virtDetPoint,
                                                sphMirrors.getRoCs(),
                                                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; }
+                 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();
@@ -584,6 +592,43 @@ SIMDQuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segment
   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
diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.h
index 3071906125c..ba850452f88 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.h
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.h
@@ -134,6 +134,17 @@ namespace Rich
               m_CoCs   = SIMD::Point<FP>(X,Y,Z);
             }
           }
+        public:
+          ///< send to std::ostream
+          inline friend std::ostream& operator << ( std::ostream& s,
+                                                    const MirrorData& data )
+          {
+            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; }
@@ -156,39 +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
-          // ( 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 );
-          }
-        }
+        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,
-- 
GitLab


From 034b0be3d2fff228a3a86a1c477c6d4a36f9b67b Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Tue, 5 Dec 2017 12:12:57 +0000
Subject: [PATCH 21/29] Tune the #photons reserve size prediction a bit better.

---
 .../src/RichSIMDQuarticPhotonReco.cpp                    | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
index f7be20440c4..e37d087f739 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 );
   
@@ -587,6 +587,13 @@ 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;
-- 
GitLab


From 41997eb438c423d337df6fdd80f13cbd98045b01 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Tue, 5 Dec 2017 23:26:31 +0000
Subject: [PATCH 22/29] Improve the default padding info to assign better
 Rich1/Rich2 z values.

---
 .../src/RichRecSIMDPixels.cpp                 | 37 +++++++++++++------
 1 file changed, 26 insertions(+), 11 deletions(-)

diff --git a/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp b/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp
index aaae787467c..c17ffcf4737 100644
--- a/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp
+++ b/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp
@@ -32,9 +32,6 @@ void SIMDPixelSummariesCreator::fill( const Rich::PDPixelCluster::Vector& cluste
   // NB : Not zero to help precision problems downstream in photon reco.
   const SIMDPixel::FP defZ = 5000.0;
 
-  // Default SmartID to fill
-  const LHCb::RichSmartID defID();
-
   // Working Vc values
   SIMDPixel::SIMDFP gX(0), gY(0), gZ(defZ); // global positions
   SIMDPixel::SIMDFP lX(0), lY(0), lZ(defZ); // local positions
@@ -57,12 +54,22 @@ void SIMDPixelSummariesCreator::fill( const Rich::PDPixelCluster::Vector& cluste
                               effArea, scClusIn, mask ); 
       // reset
       index = 0;
-      effArea = SIMDPixel::SIMDFP::Zero();
-      gX = gY = lX = lY = SIMDPixel::SIMDFP::Zero();
-      gZ = lZ = SIMDPixel::SIMDFP(defZ);
-      scClusIn = -1;
-      mask = SIMDPixel::Mask{false};
-      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;
+        gZ[i] = lZ[i] = ( Rich::Rich1 == m_lastRich ? 1600 : 10500 );
+        scClusIn[i] = -1;
+        mask[i] = false;
+        smartIDs[i] = LHCb::RichSmartID();
+      }
     };
 
   // Loop over all pixels
@@ -89,6 +96,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
@@ -127,8 +136,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
   
-- 
GitLab


From b63e14248a4cb606f9d2bb9385fa3eb8a7aa9588 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Wed, 6 Dec 2017 12:17:04 +0000
Subject: [PATCH 23/29] Small cleanup to SIMD photon reco algorithm

---
 .../src/RichSIMDQuarticPhotonReco.cpp                    | 9 ++++-----
 .../src/RichSIMDQuarticPhotonReco.h                      | 8 ++++----
 2 files changed, 8 insertions(+), 9 deletions(-)

diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
index e37d087f739..3615d7f6d4d 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.cpp
@@ -530,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;
         // --------------------------------------------------------------------------------------
 
         // --------------------------------------------------------------------------------------
@@ -651,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 ba850452f88..337e43117cc 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.h
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichSIMDQuarticPhotonReco.h
@@ -108,8 +108,6 @@ namespace Rich
             // Have the mirrors changed ?
             if ( UNLIKELY( mirrors != cache_mirrors ) )
             {
-              // update the cache
-              cache_mirrors = mirrors;
               SIMDFP A,B,C,D; // Plane params
               SIMDFP X,Y,Z;   // CoC params
               for ( std::size_t i = 0; i < SIMDFP::Size; ++i )
@@ -130,8 +128,10 @@ namespace Rich
                 Y[i] = n.Y();
                 Z[i] = n.Z();
               }
-              m_NormPs = SIMD::Plane<FP>(A,B,C,D);
-              m_CoCs   = SIMD::Point<FP>(X,Y,Z);
+              // update the cache
+              cache_mirrors = mirrors;
+              m_NormPs      = SIMD::Plane<FP>(A,B,C,D);
+              m_CoCs        = SIMD::Point<FP>(X,Y,Z);
             }
           }
         public:
-- 
GitLab


From 067b7613e72f097b008410ed93e4d826289320d7 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Wed, 6 Dec 2017 13:09:24 +0000
Subject: [PATCH 24/29] update comments

---
 .../RichFutureRecEvent/src/RichRecSIMDPixels.cpp | 16 +++++++---------
 1 file changed, 7 insertions(+), 9 deletions(-)

diff --git a/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp b/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp
index c17ffcf4737..fcaee5a16c4 100644
--- a/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp
+++ b/Rich/RichFutureRecEvent/src/RichRecSIMDPixels.cpp
@@ -28,16 +28,12 @@ void SIMDPixelSummariesCreator::fill( const Rich::PDPixelCluster::Vector& cluste
   // last side
   Rich::Side         m_lastSide { Rich::InvalidSide     };
 
-  // Default z position for invalid entries. 
-  // NB : Not zero to help precision problems downstream in photon reco.
-  const SIMDPixel::FP defZ = 5000.0;
-
   // Working Vc values
-  SIMDPixel::SIMDFP gX(0), gY(0), gZ(defZ); // global positions
-  SIMDPixel::SIMDFP lX(0), lY(0), lZ(defZ); // local positions
-  SIMDPixel::SIMDFP effArea(0);             // effective area
-  SIMDPixel::ScIndex scClusIn(-1);          // indices to original clusters
-  SIMDPixel::Mask mask { false };           // selection mask
+  SIMDPixel::SIMDFP gX(0), gY(0), gZ(0); // global positions
+  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;
@@ -65,6 +61,8 @@ void SIMDPixelSummariesCreator::fill( const Rich::PDPixelCluster::Vector& cluste
         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;
-- 
GitLab


From 100e496af5537a66ffc525e9fe9d8f7985cc5924 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Thu, 7 Dec 2017 07:35:15 +0000
Subject: [PATCH 25/29] update example options

---
 Rich/RichFutureRecSys/examples/RichFuture.py | 27 ++++++++++----------
 1 file changed, 14 insertions(+), 13 deletions(-)

diff --git a/Rich/RichFutureRecSys/examples/RichFuture.py b/Rich/RichFutureRecSys/examples/RichFuture.py
index 7cd50ef8184..1136fa64155 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
@@ -130,7 +130,7 @@ ApplicationMgr( TopAlg = [all],
                 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" ]
-- 
GitLab


From 29adeba6981d26146c221dc7f2874b19eaaafe81 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Thu, 7 Dec 2017 07:35:48 +0000
Subject: [PATCH 26/29] minor type changes

---
 .../src/RichGeomEffCKMassRings.cpp                            | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.cpp
index a842767017f..2282d4ffc53 100644
--- a/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.cpp
+++ b/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.cpp
@@ -69,13 +69,13 @@ GeomEffCKMassRings::operator()( const LHCb::RichTrackSegment::Vector& segments,
       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] )
-- 
GitLab


From 2e9d7e7c643eeaeb9dd9e5105b6bf1d9a6c7a575 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Sun, 10 Dec 2017 09:23:22 +0000
Subject: [PATCH 27/29] Use approx log(exp(x)-1) in scalar version

---
 .../src/RichGlobalPIDLikelihoodMinimiser.cpp             | 3 ---
 .../src/RichGlobalPIDLikelihoodMinimiser.h               | 9 +++++----
 2 files changed, 5 insertions(+), 7 deletions(-)

diff --git a/Rich/RichFutureGlobalPID/src/RichGlobalPIDLikelihoodMinimiser.cpp b/Rich/RichFutureGlobalPID/src/RichGlobalPIDLikelihoodMinimiser.cpp
index 3d5e2ea496b..49f6462521c 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 24db1009a3e..32ed22f7582 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
-- 
GitLab


From 875eefb4649dea7105aae728074885981a019a80 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Tue, 12 Dec 2017 08:36:15 +0000
Subject: [PATCH 28/29] update tests

---
 Rich/RichFutureRecSys/examples/RichFuture.py |  6 +++---
 Rich/RichRecTests/src/PhotonReco/main.icpp   |  2 +-
 Rich/RichRecTests/src/RayTracing/main.icpp   | 16 +++++++++-------
 3 files changed, 13 insertions(+), 11 deletions(-)

diff --git a/Rich/RichFutureRecSys/examples/RichFuture.py b/Rich/RichFutureRecSys/examples/RichFuture.py
index 1136fa64155..e1c09378588 100644
--- a/Rich/RichFutureRecSys/examples/RichFuture.py
+++ b/Rich/RichFutureRecSys/examples/RichFuture.py
@@ -118,13 +118,13 @@ all = GaudiSequencer( "All", Members = [ fetcher, odinDecode, richDecode,
                                          tkUnpack, tkFilt, RichRec ] )
 
 # Uncomment to enable monitoring
-all.Members += [ RichMoni ]
-all.Members += [ mcPUnpack, mcRiSumUnP, RichCheck ]
+#all.Members += [ RichMoni ]
+#all.Members += [ mcPUnpack, mcRiSumUnP, RichCheck ]
 
 all.MeasureTime = True
 
 ApplicationMgr( TopAlg = [all],
-                EvtMax = 5000,      # 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 )
diff --git a/Rich/RichRecTests/src/PhotonReco/main.icpp b/Rich/RichRecTests/src/PhotonReco/main.icpp
index a13821b4941..670956f1848 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 83568996528..0d3729b9b31 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));
-- 
GitLab


From bbacaabbc8e846900457e697b67547ced66c97c3 Mon Sep 17 00:00:00 2001
From: Chris Jones <jonesc@hep.phy.cam.ac.uk>
Date: Tue, 12 Dec 2017 08:55:42 +0000
Subject: [PATCH 29/29] Some explicit float type fixes for clang

---
 .../src/RichCommonQuarticPhotonReco.h                     | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h
index 49c5b167d3d..80fc56cc92d 100644
--- a/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h
+++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichCommonQuarticPhotonReco.h
@@ -162,14 +162,14 @@ namespace Rich
 
         /// Minimum active segment fraction in each radiator
         Gaudi::Property< RadiatorArray<float> > m_minActiveFrac
-        { this, "MinActiveFraction", { 0.2, 0.2, 0.2 } };
+        { this, "MinActiveFraction", { 0.2f, 0.2f, 0.2f } };
 
         /// Minimum tolerance for mirror reflection point during iterations
         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
-- 
GitLab