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