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