From 3ba3a9de8462a4f34c323e8df143ba0a458c2a05 Mon Sep 17 00:00:00 2001 From: YiHui Liu Date: Mon, 4 Jul 2022 23:52:05 +0800 Subject: [PATCH] Initialization --- .clang-format | 208 ++++++++++++++++++++++++++++++++++ .gitignore | 7 ++ .vscode/c_cpp_properties.json | 23 ++++ .vscode/settings.json | 0 .vscode/tasks.json | 30 +++++ CMakeLists.txt | 19 ++++ Q3D.code-workspace | 106 +++++++++++++++++ Q3D.exe | Bin 0 -> 78848 bytes include/GaussNewton.h | 89 +++++++++++++++ include/LevenbergMarquardt.h | 121 ++++++++++++++++++++ include/clip.h | 92 +++++++++++++++ include/getADC.h | 192 +++++++++++++++++++++++++++++++ include/utils.h | 20 ++++ main.cpp | 6 + process.cpp | 9 ++ readme.md | 0 16 files changed, 922 insertions(+) create mode 100644 .clang-format create mode 100644 .gitignore create mode 100644 .vscode/c_cpp_properties.json create mode 100644 .vscode/settings.json create mode 100644 .vscode/tasks.json create mode 100644 CMakeLists.txt create mode 100644 Q3D.code-workspace create mode 100644 Q3D.exe create mode 100644 include/GaussNewton.h create mode 100644 include/LevenbergMarquardt.h create mode 100644 include/clip.h create mode 100644 include/getADC.h create mode 100644 include/utils.h create mode 100644 main.cpp create mode 100644 process.cpp create mode 100644 readme.md diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..0bddc96 --- /dev/null +++ b/.clang-format @@ -0,0 +1,208 @@ +--- +Language: Cpp +# BasedOnStyle: Google +AccessModifierOffset: -4 +AlignAfterOpenBracket: Align +AlignArrayOfStructures: None +AlignConsecutiveMacros: None +AlignConsecutiveAssignments: None +AlignConsecutiveBitFields: None +AlignConsecutiveDeclarations: None +AlignEscapedNewlines: Left +AlignOperands: Align +AlignTrailingComments: true +AllowAllArgumentsOnNextLine: true +AllowAllConstructorInitializersOnNextLine: true +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortEnumsOnASingleLine: true +AllowShortBlocksOnASingleLine: Never +AllowShortCaseLabelsOnASingleLine: false +AllowShortFunctionsOnASingleLine: All +AllowShortLambdasOnASingleLine: All +AllowShortIfStatementsOnASingleLine: WithoutElse +AllowShortLoopsOnASingleLine: true +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: true +AlwaysBreakTemplateDeclarations: Yes +AttributeMacros: + - __capability +BinPackArguments: true +BinPackParameters: true +BraceWrapping: + AfterCaseLabel: false + AfterClass: false + AfterControlStatement: Never + AfterEnum: false + AfterFunction: false + AfterNamespace: false + AfterObjCDeclaration: false + AfterStruct: false + AfterUnion: false + AfterExternBlock: false + BeforeCatch: false + BeforeElse: false + BeforeLambdaBody: false + BeforeWhile: false + IndentBraces: false + SplitEmptyFunction: true + SplitEmptyRecord: true + SplitEmptyNamespace: true +BreakBeforeBinaryOperators: None +BreakBeforeConceptDeclarations: true +BreakBeforeBraces: Attach +BreakBeforeInheritanceComma: false +BreakInheritanceList: BeforeColon +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: false +BreakConstructorInitializers: BeforeColon +BreakAfterJavaFieldAnnotations: false +BreakStringLiterals: true +ColumnLimit: 100 +CommentPragmas: '^ IWYU pragma:' +CompactNamespaces: false +ConstructorInitializerAllOnOneLineOrOnePerLine: true +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: true +DeriveLineEnding: true +DerivePointerAlignment: true +DisableFormat: false +EmptyLineAfterAccessModifier: Never +EmptyLineBeforeAccessModifier: LogicalBlock +ExperimentalAutoDetectBinPacking: false +FixNamespaceComments: true +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IfMacros: + - KJ_IF_MAYBE +IncludeBlocks: Regroup +IncludeCategories: + - Regex: '^[<"].*\.h[>"]' + Priority: 0 + SortPriority: 0 + CaseSensitive: false + - Regex: '^[<"].*\.hh[>"]' + Priority: 1 + SortPriority: 1 + CaseSensitive: false + - Regex: '.*' + Priority: 2 + SortPriority: 2 + CaseSensitive: false +IncludeIsMainRegex: '([-_](test|unittest))?$' +IncludeIsMainSourceRegex: '' +IndentAccessModifiers: false +IndentCaseLabels: true +IndentCaseBlocks: false +IndentGotoLabels: true +IndentPPDirectives: None +IndentExternBlock: AfterExternBlock +IndentRequires: false +IndentWidth: 4 +IndentWrappedFunctionNames: false +InsertTrailingCommas: None +JavaScriptQuotes: Leave +JavaScriptWrapImports: true +KeepEmptyLinesAtTheStartOfBlocks: false +LambdaBodyIndentation: Signature +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: None +ObjCBinPackProtocolList: Never +ObjCBlockIndentWidth: 2 +ObjCBreakBeforeNestedBlockParam: true +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: true +PenaltyBreakAssignment: 2 +PenaltyBreakBeforeFirstCallParameter: 1 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakString: 1000 +PenaltyBreakTemplateDeclaration: 10 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 200 +PenaltyIndentedWhitespace: 0 +PointerAlignment: Left +PPIndentWidth: -1 +RawStringFormats: + - Language: Cpp + Delimiters: + - cc + - CC + - cpp + - Cpp + - CPP + - 'c++' + - 'C++' + CanonicalDelimiter: '' + BasedOnStyle: google + - Language: TextProto + Delimiters: + - pb + - PB + - proto + - PROTO + EnclosingFunctions: + - EqualsProto + - EquivToProto + - PARSE_PARTIAL_TEXT_PROTO + - PARSE_TEST_PROTO + - PARSE_TEXT_PROTO + - ParseTextOrDie + - ParseTextProtoOrDie + - ParseTestProto + - ParsePartialTestProto + CanonicalDelimiter: pb + BasedOnStyle: google +ReferenceAlignment: Pointer +ReflowComments: true +ShortNamespaceLines: 1 +SortIncludes: CaseSensitive +SortJavaStaticImport: Before +SortUsingDeclarations: true +SpaceAfterCStyleCast: false +SpaceAfterLogicalNot: false +SpaceAfterTemplateKeyword: true +SpaceBeforeAssignmentOperators: true +SpaceBeforeCaseColon: false +SpaceBeforeCpp11BracedList: false +SpaceBeforeCtorInitializerColon: true +SpaceBeforeInheritanceColon: true +SpaceBeforeParens: ControlStatements +SpaceAroundPointerQualifiers: Default +SpaceBeforeRangeBasedForLoopColon: true +SpaceInEmptyBlock: false +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 2 +SpacesInAngles: Never +SpacesInConditionalStatement: false +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInLineCommentPrefix: + Minimum: 1 + Maximum: -1 +SpacesInParentheses: false +SpacesInSquareBrackets: false +SpaceBeforeSquareBrackets: false +BitFieldColonSpacing: Both +Standard: Auto +StatementAttributeLikeMacros: + - Q_EMIT +StatementMacros: + - Q_UNUSED + - QT_REQUIRE_VERSION +TabWidth: 8 +UseCRLF: false +UseTab: Never +WhitespaceSensitiveMacros: + - STRINGIZE + - PP_STRINGIZE + - BOOST_PP_STRINGIZE + - NS_SWIFT_NAME + - CF_SWIFT_NAME +... + diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..c60b80c --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +# data file +*.root +*.log + +# build cache +build/ +2016Q3D/ diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json new file mode 100644 index 0000000..ca9e844 --- /dev/null +++ b/.vscode/c_cpp_properties.json @@ -0,0 +1,23 @@ +{ + "configurations": [ + { + "name": "Win32", + "includePath": [ + "${default}", + "${workspaceFolder}/include" + ], + "defines": [ + "_DEBUG", + "UNICODE", + "_UNICODE" + ], + "windowsSdkVersion": "10.0.19041.0", + "cStandard": "c17", + "cppStandard": "c++17", + "intelliSenseMode": "windows-msvc-x64", + "configurationProvider": "ms-vscode.cmake-tools", + "compilerPath": "D:/Microsoft/VisualStudio/2022/Community/VC/Tools/MSVC/14.31.31103/bin/Hostx64/x64/cl.exe" + } + ], + "version": 4 +} \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..e69de29 diff --git a/.vscode/tasks.json b/.vscode/tasks.json new file mode 100644 index 0000000..57f0f7e --- /dev/null +++ b/.vscode/tasks.json @@ -0,0 +1,30 @@ +{ + "tasks": [ + { + "type": "cppbuild", + "label": "C/C++: g++.exe 生成活动文件", + "command": "D:\\Scoop\\apps\\gcc\\current\\bin\\g++.exe", + "args": [ + "-fdiagnostics-color=always", + "-g", + "${file}", + "-o", + "${fileDirname}\\${fileBasenameNoExtension}.exe", + "-I", + "${fileDirname}\\include\\", + ], + "options": { + "cwd": "${fileDirname}" + }, + "problemMatcher": [ + "$gcc" + ], + "group": { + "kind": "build", + "isDefault": true + }, + "detail": "调试器生成的任务。" + } + ], + "version": "2.0.0" +} \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..fb16bf4 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,19 @@ +cmake_minimum_required(VERSION 3.16) + +project(Q3D) +set(Eigen3_DIR D:/Microsoft/vcpkg/installed/x64-windows/share/eigen3) +SET(CMAKE_TOOLCHAIN_FILE D:/Microsoft/vcpkg/scripts/buildsystems/vcpkg.cmake) + +find_package(ROOT REQUIRED) +find_package(Eigen3 CONFIG REQUIRED) + +include(${ROOT_USE_FILE}) +include_directories(${PROJECT_SOURCE_DIR}/include) + +file(GLOB sources ${PROJECT_SOURCE_DIR}/src/*.cpp) +file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.h) + +add_executable(Q3D main.cpp ${sources} ${headers}) +target_link_libraries(Q3D PRIVATE ${ROOT_LIBRARIES} Eigen3::Eigen) + +install(TARGETS Q3D DESTINATION F:/NuclearAstroPhy/Q3D-Calibration) diff --git a/Q3D.code-workspace b/Q3D.code-workspace new file mode 100644 index 0000000..ac2245a --- /dev/null +++ b/Q3D.code-workspace @@ -0,0 +1,106 @@ +{ + "folders": [ + { + "path": "." + } + ], + "settings": { + "files.associations": { + "Offlinecal201609副本.C": "cpp", + "Offline201609.C": "cpp", + "Offlinecal201609.C": "cpp", + "Trans.C": "cpp", + "Trans1.C": "cpp", + "Trans2.C": "cpp", + "Transcal.C": "cpp", + "Offlinecal201609_Li.C": "cpp", + "Offlinecal201609_d.C": "cpp", + "hadd1.C": "cpp", + "hadd.C": "cpp", + "Filter5.C": "cpp", + "Filter4.C": "cpp", + "Filter3.C": "cpp", + "Filter2.C": "cpp", + "Filter.C": "cpp", + "Er.C": "cpp", + "xstring": "cpp", + "array": "cpp", + "chrono": "cpp", + "compare": "cpp", + "format": "cpp", + "memory": "cpp", + "ratio": "cpp", + "tuple": "cpp", + "type_traits": "cpp", + "utility": "cpp", + "xmemory": "cpp", + "xtr1common": "cpp", + "xutility": "cpp", + "span": "cpp", + "vector": "cpp", + "cmath": "cpp", + "iostream": "cpp", + "core": "cpp", + "dense": "cpp", + "algorithm": "cpp", + "atomic": "cpp", + "bit": "cpp", + "cctype": "cpp", + "charconv": "cpp", + "clocale": "cpp", + "complex": "cpp", + "concepts": "cpp", + "condition_variable": "cpp", + "cstdarg": "cpp", + "cstddef": "cpp", + "cstdint": "cpp", + "cstdio": "cpp", + "cstdlib": "cpp", + "cstring": "cpp", + "ctime": "cpp", + "cwchar": "cpp", + "deque": "cpp", + "exception": "cpp", + "forward_list": "cpp", + "functional": "cpp", + "initializer_list": "cpp", + "iomanip": "cpp", + "ios": "cpp", + "iosfwd": "cpp", + "istream": "cpp", + "iterator": "cpp", + "limits": "cpp", + "list": "cpp", + "locale": "cpp", + "map": "cpp", + "mutex": "cpp", + "new": "cpp", + "optional": "cpp", + "ostream": "cpp", + "queue": "cpp", + "set": "cpp", + "sstream": "cpp", + "stdexcept": "cpp", + "stop_token": "cpp", + "streambuf": "cpp", + "string": "cpp", + "system_error": "cpp", + "thread": "cpp", + "typeinfo": "cpp", + "unordered_map": "cpp", + "unordered_set": "cpp", + "xfacet": "cpp", + "xhash": "cpp", + "xiosbase": "cpp", + "xlocale": "cpp", + "xlocbuf": "cpp", + "xlocinfo": "cpp", + "xlocmes": "cpp", + "xlocmon": "cpp", + "xlocnum": "cpp", + "xloctime": "cpp", + "xstddef": "cpp", + "xtree": "cpp" + } + } +} \ No newline at end of file diff --git a/Q3D.exe b/Q3D.exe new file mode 100644 index 0000000000000000000000000000000000000000..84a4999ff237f503c6fd8d1b6ddae41f94bad0c6 GIT binary patch literal 78848 zcmeFa4SZGAmG^&>8z2I5ZYoh|sm3<8apJVG*cN(en_!~nz&RQfYEng0Ul7N)s7b{! zfzXha+ry!>Pdmk#_9-*d{%dEP(V538c5IsvUIHitR1EkQ-*1dB5i0}~?*IGS=iEy` zblQ1lKF|O8{D%*@XP>?H>)LCtz4qE`uj74vWh@?x#R~Y#WMZ*3JoT^A{C@b)FyINN z{ANPzzT+M_bxof8$fG$Py~_G&{#VTNC;8#~-{gPVyx*CB z8qY6IxikM%o_{~#&iqoIpS{S&uil0^RLbOh1XuWKy_q$R9h2^UG?qJ zv4scMUXhI(ihUwBE`P!)u|338hq13IoG{_aJUtZ7kHtotsKd`#8&DyUfALtul5D0Z zFR@`E*@x0aF$%`c`dnVDKovCP#rpNE&5JdS5hVK8m={|>F!8ax*pzJCh-Vw|WxV%J z&WnW&jev=~*bwdgpnsofx_V&~z|RfwqXwzJQ9lToe+HaqjcbO(4Bl`lOwjY^UlA`Aayz&xHRLdmWukCwd*5nu;8MqkX5dZ@-)Dw7%EISUG;LW8Jeo z{ps(=Vk_<7wn>hEBf0N0#oDQp?bKz3cE{^YAD_A@Q~Z1^FUHHF<8A+ubPeQYr$!IP zO)4@4x0QWZ=3k7=TzGWm)aBFWZ#)-r@%(cN9e=c)sxB>Y`6=`K^Ahg#a|@eCyVK84 zH2v8JbmpHxgX3?fGrhLIEq(DyVddOH+qZ@(u-LYfQC#+N8QMONBD^bvA6&38V+q?7ab$QigmtQuoJ3Vm+!1M%u-#^|A zK~rf$&6!m7HO9Nd_HSWyqft-ZT-bE7?H`vZu4VnOJN9Zswx<8rC|*)yG$N2ZO7TkjaK_IlK2}=vtqHf z)-jo$_MV|kXSKhzy1eJX5z8-dok)N&JmHCTBJQSr2Vi)}k;~Z&h+v%iq{nZn3Z2 zX(u;YFe!+KLdNk7?|J|9X zg6}F(RjTmzXKepc#$ChE!qEFoP(u5?1K)|EIYonW@@MZXr_2HrL0Wv-x-UkUyi>m zQ~cz+L^=K#wSSiHP9K531nAqK;atE7R%j|!U3z&mj@enH@nrs7uR(1R=cHOQi7Flc zEgDu@qJ}ws+RUb1HB&dy6!`Cq^1Rr;=A|#CDYienMk7wwo5uq>6;S%i1~AU`w{C!q zz+KjhJvshPQ-kY&-$p{Q{Vi3Q9cp2%?O$IC>gm1cP7SU!Bp1=!^=cXfBAS~UQC+R) zL*4^cBHt>aY%4JZ~8ff&7-{OR@3Lj3r`hepJJQ#(&WX5 z`TO@E2<3{{q63Sl^kXr4>-txfR(KsRS}Q(InmtzH_=t(qqw zWyh-7;`m**bx+qG80K0hn3-Yq8Id##@seu~$2$Jnky-MNoyBSIDsbA@G6fH}$4zdv zmfr*>+~gyzpEr;lN3@L=?VxczfLJY}R>50T#tt&IYZsV4COmx~ZCyvi(L{c6*^c{Z_;f zd-@%TCYM#W>ErN+rb)T+V(QuZGInyv`rZV?3Z8(3ce?&IS3GM=x|_DSzW0FZ&%4L< z=l+Z)eb2PQ_V0jUwX#03bk{s<`o%CwJSAuaWtP@f1g*r&uM@8kQd+QrkZb?y4UN-SS$XD zsju}s#7Ho4h&->jzC%~@W(8By8n-f-OLL)h@s0sRRT_VDAqBh%!}!P^tNw)F6Bl!n z+pHBQ5IQfI_>Pe8wC^vlR%{^T`NK}vE2BLBWvA=KQL}>r3#gGyl)4bM`QDIur33m19jV-_+w6do86wyYktu6o8JU^uV_A% zx&95;xPB`mclBZ?xxe)b=LJ*Nt1^VcjkD+kV*~TppnkaJ+eX&m#CJIHEs{N&(|f9{ zqKSX7GYo$;)J}4rwfwu{Djej_P4kMO4WCcH@!`8)5I`NHeFk`K&h0Ho7 zW|ch(3pZF^{xVJS{C7hvKQAoPxOr-*6Mux_ehK5~A?j;{e+&wHJpWaL^zI1htY4e$ zeY@=VCI^nQ1E6h z+DeQIla9$nDZ!IOb}Gvg*=Q&Jo)dr848{KxK8&=WVSsnHj2WB@Aug(gEFV{UW(6~z z*3(1UcjG(NIXcl;4=p-s z?jXX$utv_Ee^uA~A*bu*QHYTm-0a*zU^jEeWlHUP%Jb6C6+ zlbtWmn>jHM%@^iMJ#!|SOI=WL{rG6UJmJK*=H?4TOOI`T9&F==Yd&b6zy$xjd6Gr! zp0hPy5U|Cj4x?Rj=C6*LGYs-JJ6^MI1c!84Op?Q{OgvYY)E$R|P2->NNg`io2SdQM>tn@S{^Nh;p|jGKC>G(oV&f7hFSS;AT| zix+sdVSUmH)agCtHP|F9AuLAh`0pY;SSt?EW5<67$$d^~Ax&dNjXQo1!@o_t?*&>< z8m)w@P)q4THyGc;V%z5U1*NvX*(&nxa@mB{S&3@);HjEr%oJ}KgRsiXD@8AKQ&*Kv zqx~@4BCkhfRghqfOUad*zqXRr0_$f&Y4W?JwJuR5N{{G(>`(e$1u6*)v8hMoyts4l z^hiJ8B1)+U#vE)98BElL_~XEq9VaE?L<&RRZJTLs^1a55X+_<)yO zl(3e+Y#2v~CEE`%=lUV0T|Y!KWTMg6CgzfP-5_dvM0W^~6KFzhS~@zqF{*B_c!bCP zV4X;>7WrD@{HdQbv^JuU^zXWk7BG$cWwt_xMKCl9SKEHUX?2GhE7`^=u zYGbMO*WYuVKe2bFUr;oC)VF|H-`^hu!tuD92y#ly9*2F)z?poSSjtT}IAw5gd{|i%&#uSN^!k9}n{ne$sC5 zg!ncc331CiouD%As3~rbof_vQU$RzQg}&7FCrsVpmA{122J*5Te;rlFKd92yPuR#f zEmMwFwVj$U@DWN>13svBSQMl`Vc=D;#U|g?JozdU*$5RhEV-ZdbusZPtmPZB7-d6S zYv1YmPk>YAPW~DIF{f+)Xs7+TVK>>MLSey+zCtOUzdhUrvWx22(R7-TB;f>_${H3n zWKMafodNSUv0eRfx4hR&okUAq|6~Wj%-FBu)1_>(W~ zTGJ=8@3gWbBtu%cHrvE}w;ZI68@N~PPad#3eo258f34?lap-}wWLR4C#SqLJgD11O zZJMkoTD>??^;ZjqY^H4TTD?G`iDGjg~5aIeTeVA0X|MbQfk%DE$IQ3F*}vnJ-^1 zCYZ50{>OW;l;HTGj{=xI70d?5zo+oWNHp+GfLckye}CNko}n7W61BAfMX}caVn&Xxg3qmfr6H~b});r-9Wxn=|MrLQ7Z2>yXnLG1F2$C zW6bZlX#RN7aflw1?DmZ%wsy8p(w0S&vikR!o!LwlpoPRQenZQRsnswj*MDkmq*q7e zU53cJP-~iBS%Kui{><@_LieKUtnmD6P(DT;#-N)SN^rKrYoP%>hFQ)MMCrop!r5*TCq&#-+-oL|9OwgE|M3KGax$3J%xV$`K* zX03bdNmkL^d#PkKDrb_OTVki0S4v>)mHHEvQXh?=F>{x(E9swa-MYI_*xyc;PTQ}d zwKZvWBW(YwIVwPTN<^Ea$jT6^b0v4Hdu z;sph5S+tSvU>^Lm^$XN%cn11TiRZtr5BB8635uT{3$v~8WVAIDi`=$k-d zx9y)@MolO-R2gODkx6;U18SIpLaKCYn(8fNuJv9;DGRLpooh{hohe zp_Q1d!H4Ln3A$O#KqL$eXdWQe6f8PV=NQ4o?|S9CnRYNn9Z*ylj+mt%dOk!ZH$Eg4 zl9d26b-RcdA6L&6Ck2?@{%T(HV}y3Q$&HJO?XDL_kz}%TEhj{2cF{)&N`uvB2x>JA zHWag2(zs~pR&CLPnT1YyFT_aPs}PcD7KF4i?At_C3!LBri0RdQ_F4%L%2xbG%8|vVX{*|AoAygNOKCI%iDwxR zoSWD_D-ox_H6;>X&di-LonrfcnUK=yJ6N5uGixo0HL<<)N(Sp#5)P7?MVRx!{Mc=u zXI--Fz_xUPIM#k}2rU)=EW7==Av721IsZo5^5=}j3hRK0AN)wB__e3hq#o(?mdbnt z)@sT{NOPw*m$+8-TQuonEIA@csMkeX>B-Q>iyN&Lyc?Vq9)vAHC8~iyBlckfvE@M;p@EQfO8Z6z;+3JnVuOcJ8D+UQj~omAD34%~T1 zG%@c`)>`Rx!q+^Yc_X0cG`kACfc2JBbtk!{)xq^Aw05O|+ zSyN^d^cdI8qxLVE@pt^%#E6NR?6to0F=_s?y0JfU9zOg4dCH%2HT)bzS&oMPq)hQu zyQG&4$G)4fmtH-Auof!bQZkA<3l8>?XgdeEMAMfgt%)p?f-)aPdV)!LW#Ult?-yA@Lq0U*hcr|M4Fy+M zW>tNVW63|fj~qDTI zEKebKdm)b^6LFj(s<$T;q4->qNs1A`x=d*Rq8-5!DcRfdFVxs7_~9}sOj`?pY4&-q z6`D4krnQ0wsWtv_G|2NWfa;G^yR79esl{xa6K?SJd~}~zC zk?kA5L#^7LEhDTTi?U}>CKVo*;FY;G!kgOXZTj# zTZo2bJAKN&BbYFLhu>jH)!`sI$wylrR5cs~j4Ue}!dDNs#q1|IEK1c;dLtz)EIGWo z1b5jw1YGjkN>&By^u`{~e~u%2V|K_X36f4#!px2J ziCB@qg^|5l+oZ^{5;dHX)~w=$wE;!fT~^W6YZ!oEC?|}km*dxdtLR6od9{j+)J2XG ztEi(>F$u5DDhiCul&D!rY=0QLn%E%|tU;gF5ucA|fJm`XZ<>TH>kdV%Wnm@2(~ zs4y_14QH{mbFXONtWq|FUjXujfjqiy2c{`@;4Mq9m@HKM=w9k^Q*-Zi{hx+%^}IER zi)s_gN)LJ~)vk4suE59H!_Hl4L>NvZ{dvo%`AqhmXA%F~K4szvbJ17MOT8yb%v}PJ3QWyPQ zBgHOLU;Md=^b<=$ZW+rIpP?2F=wR1wUsMSFTwxA!zaqi`I}NNh5YAqEV@DggX_qKw z0wX5dBkbztF;QyDkSWNz0I*15YKDHSAs$rEPJ|L3ROm%ixWi8UNQE)9tIZ>%p)CCv z^3Yg`JV?`Tf+ER|p)j_ASVKQWKN8J=a??1vZV1Or5>CR7pp2%Ij-U+N|B)yIau;Qo zK~AC4O?iE;Z^)GOsUJc?!pjm$8YLuGB;H3X;3G#Z!zxqMQkkU|juHz$gj`nA_^d(A zG3aH*w!`#NWz~HV(kc4lz&M5|q!@I$7$VS7uPyz96RF#0*OaQcizn=+h6wFf+E~q5_jt?fxTlu1+5Y$*J2n1*fpNHSOt5>P z(7@Zf@;@j6q(tFTN*+}LNG#|P!*&@qHKnO1G3gg4!591 z_f}?#Pu>c?(b9{IAviq^4@_7^%_Sq~tOn5yG%@EiRdXBjhuO7?>P%!{_NVE8d^?jt zGBHT~J!Z8l2VAm56MmYRDz=ANnAp4X8yc;Di=pABk4Z*jk>QYAr}J>^?|spbW7pvE zE_*IPl?0+3kz^uCvR7RKZ(91gc~J03mNG5fM-;uB7lWJdAO)3n*Iko`;|6UkNt=F%-|SxRnCr8yPz87O zg==$yVp*H}JpV0>kzMTBWr6ik*cwYOOK~}TlWr#%25afbykMHe3leaVII7TS7UdrZ z7d}i?lSlZYuM1(V%OW~h$z)ya`eWd@%$SRT&%GBO+G$qixg5^a+yy1BL4yjLMH@t6 zYFg&`@Dw@Z(Ix(*g%}WJ?3#rbnlz~$tDQYRQja0DH zxs9Z%SuGCQMk)aP@K~)EL=O?WDvXU0dsi6Cl;{u(B`s?1iCcj5gpdj#;;0nqRZ?V- z<^@Tm`6PO)#cwn}bz{WK#taOAVfBRN44W_vyXXSgn7f*(WHXhf(OEOs(pS8i$sC>g zAIXglO^N6bZT4s^bWN-Ks9C$1?>Pti}Ku?dGe92*F>K@<<^M9mzgPc&IgtU)9hF&18b`nH$fKNKYzgNU;F zX)+Za9*osyFo*?3P%Ui+M>Qmws|7F|tJNr1G(Kj?Amf0wB|Bsr82f?Wsf%!|FGQ1u zF=6#ian6cRSM?xSVRXC�FkZ`qW+8SxDZKeyT=Ed$lepeV=X+iHR9zcB0FibNy3{ zT_P;1fublYvN@#f@8e(6Qo+CGh=0iW*n6ITeSZ_R;$V!;coMVxu7L7cwI_k*-A;Z7;_9OYhJm^1oi_VeP;4 zMpFS@6HoevfnW$~=(S7T@=tJl^f6qM^GAEhzUFI5sHwS0lS{kRy0QEejU70Ryj-{x z6|=3YANxCD*iOFH`pK+t!)YVHf<4h8O^;L|aSZR?YDphmSRa`+Rl-722*jxx81y-l zcxr#gH|E$k*Jz*uRkwPFUb?2|)vTAXW$+4}@E zwP0vOhjq$GhhI6m!>`b=?0g_#<^$7sZixuq%!i6&%!hd+<^v))cCZpC%uMhWrq3T> z8mypTA?(a>-PSEJZjnIoccnAq4CS3Yi9z=K%ek(5isLVuq*b?pT$F@B2D7ajNzxQh zDfB_kG6)elYY8v${0WSH!hxi9%PW1wH*#dXZ=HrgWc@lJL)IeB49%om{5K6b!O)&~J6VQ38zTpBNr(5HU@w@klkJh>sp4(&(ld zH+6scEz@<(B{IeDVuK@b&dwAM5j34|qrG$l2szPwp=^ z04r4}{T+9>zp?+U{-Wpn507v86Y4JnJ#uh)uQTrr8cU5ZXAS@E2s2HxBmC>{9`5nv zwYd?VK=8vySezdoCgnS%RCv)`#K-V3$>mO%BZZFa_Oal%3t`!A>W8JP>AfVuT8^lv zJFpL9!+(4s+!mI-qJV&lp3l$sWLI@xsn%;$nbh)7X|D<2KFLk}ytGH<_*W!kqTi`7 zY1f#zTY&2xvv8pXE>tlR7GxsNWM5u7NsN^8RuNgXw`+{`ZCxGU(wHYRGF|@Ct=H(t zkP=GNyJ_gtYo@+v?uT7Q5-&I_PX^cF4XsLZ--WAdh35q`Dv;=z5|1z_-Xt%cLG4D30T6YG<3Gu@FE5z+uH(P$ zaP7`%Cg~e&|88>u!bz@e<)rhTJ~sEyQc%H042#GPGBznE z@D|4XC#8jCQsJYW%sSC*CE8IbON%he-T;5OqtHvfVlD5H=moe$Q*N5;N5jJn#|($Q z+7%_QL<6}yCi5;~m(GO*jNLKeiR+Re7rLgd zZGSFLy3y#`;L_4%WO9<*kNQx0(%m&_b&f$#SYbg=AsyL%h~5vOC$t*JqNgM1$Z3C#y~v}6I>4vy z#?|Uh9+bpv7h%;&4^-x&tzlkQg@_ zE?CiFi5kZ^GHjk2VajHsGa8LkzQHf1lij-i)4qWtjDyzlE;Fsep*Dl!CEw(}vcWV0 zhv&ayM^oL2ujhJ+=C(Kn76OFVmMjW zs@;FYS^tygLAHG$=!51%R<9ZI*Tla+AO5I$VEE|=(m$jxjJ1J>9HrxABBXh^l$f6; zgdGoAo-SW$VT3!vNwpfnAvlXp^z8flI5e>HNlNPVc4W$q9NlC1_tM&fk)jKUB%oTg1MQ=y+HNK z4$aZie`{nw8c~vUrR2ZHacl*K==R1u)5?uirnTU6V^u;iaL2z}^6t4;-`H$*JO?!_ zIdJWbS36dR_Ki#S-`I5Z^;X9w!bV5_Ljv&ADt~{K|5lCuZ(q5lfXFItIR0I=-&gH- zR<*xiRRtGi_N_fg@y(w9XjRwnsICL!{q7~tXX3F~SNg=v+Ds<~2uohbV2_>Y4yy91 zxEyEQ^C*#3NGweg&F|PMHvP@x`~!ny{5^y5YJX2vW^IjB1gN27RqY27ElUOvaA{nA z%ZcnXQEJcS^xvyY&#z{){nMkGJMdRTEP1!8xoJUjV#$H4Z*01@xp2w;E33Tfrb)zm z)z>!HEje)O6<0Sk*Dl%LboDK#HrrD-NgQdj{*I1+nbW$qTkauhl!>bL-&&(6P||g1 zlJ$4~fvT>T$982t;=k2(D6!-f(_JE7IbrbRYM)cL5iPuLCO(J!=1RYLn%`ZX z^qBr0e;Tb94F(r!6KWOZV=`Bm{_Sc}STO2~1{xRrL-uAO5__g%p%3~^-YPiss0k(^ z1_oV6OjahyQT6@$fn(L@MfJUDphwp?px7hp+h7v@5%tv^v%aT)Mu{UjHG;t~A(-zt zmwFo`icgH$!cK*(Ba?oxYGmhW_Dd>-5UI__0^$q+h1B* zNlAOEt;7#YZ9?Xft5!HvH8d9-O)1qTl*K2k6yFv3+K_ev)MeiGSKi|uD?3wlCA?(G z1!1OLC_0kFE@ypRZhD=)DFU?WFr;#AFK)_I1tV&!-O^tpkyayo2;TQ~KS(l!!g{bL0M zr(^bGrYS; zA8?Y7Sb^?-3%USkjaByQ{2FWG>x|FcrT41)+y_8fVc*A9<$gqc+}J7F2)gTp?;c8X z;@xgg+HDK~v+|tdKU`%^l1Uv(vf%PVj(?&PT%OTxmqowNx;RsVIe>LK=7}>jiekVdLs#x&{l4e-mmZm9Yx35~o9HHz%yJ^*D?xtVZ@na&An#F~8cX z8^*U>|z(Mab(_J~X^_!%@<;ik{nRO5w5%W!G~xqy9Bz z?$pmq6jc%>?l#r*S1Z0ROx#7d{tu>>{;(;eRViJjz9EygCydJ0D{xrv(1+Jc=c5*> zd$C--jKo!?2|D?WT2m8Ec+E5+mTSVg2vB!a-HKnU_P?*Q)FoiI=^QoSqk$Ug8?CbH z(uCjS1_`S|*3k`8{0BCun|5$} zEoxT%Q&ct-wnw$9Mv;eRftC1G8)Qz4zb_eI4WW@pN`%^OC2cU()tVvdB-v}AQ^`Y;X`mS+Tb}sU zYBh2M%B<7Dr4VY!R@4b~XWM1+hExM)A2XP3Pz$o@!bUZYY?Rud4t1K=sXXIrRTDVu;J5k$doO%pa4VSSJ4-T1{k3_9PdNfvXGy8VcJBUR55IsDK8g zG93+TGm#~VWF|5TBaj)hdd9l`5T&eFYxaetq(Ne2Ld=tesl)Fs{`#$se`cRaC|ooo zG-y#89Hms1c-*rB1i=^!uIF#&5}483ByLY|_GB8fH%9{_}8&!&F-uRIP8Ix&9FdwEl5r;c^~2C88cTRV$O;Ou@az>ycicR?=aG zUOLArnz5T7cQ@yD6L(izMRN!xst=farh!;y6zoxAgK@M$)q4^uV2Jnw3UrdI-~d0D z9!T(8ieyyOv(=t`53|W@`Tt+cIlVSpv}A^IY}mv!l71s??QC% zTnQ-13I|Sb?jbwL#w9GkTHd=)B3V{I3`HvAZ3daa@!xcl4;g9A3KmJkh?I1p;@96G z+pap6@QqTHNXkG`Gx7wt15~3B%zDg5-pC_)p&au3c|qwffC%>!B@Dc7`$Hr3eNIp? z6>+Uw>%4!m>py~%j`IqUU>@i1X*a&@Y(#-AhyqU-F`yE9#0DGt<+Lt1&tHjF1ZPu; zuX7@uxB{J0;xBY?drJvFe_1AM&aBYyEEExInM7Rn70aE7OR@Sq{w^@X)2pG`Z7y6+ ze395H(Q18#Vh!j?s|;Q^@pWS5z_y+x70`D5JY_i6#J$SZz|Xafxm-`1xiAqIy7j3x zp}?yDmz8>E^FAYT1}~jJmaGm_JgBiYMZmwrH;jfc)fvXt=UC~pfW`fn^hc1CxMIZ3 z!Wv*IW`Y>@V)T<%$0O(<5&X?u9wq~(V$p3xIdHqu{ zdrYK)YAs?BF1Z@qmP6N51}*w|TFO?dqmkL};7<8-GI%T3E$OkerBGVg5IGuQ9IT1= zkjo5ZsHc4l4RQ4rv`SM-WgrMHALH*~@rFm8xI{JhIyt0B-a?=US7YP1@H1x%hj-xw zG&M!m_q!`Zfclf1XQZ(1%l`-5U2b+Kp8uND^-_Vuet@f5 z+J5vZb5LMw2QWBcgua==&9>xGG*em&}R(sO2Ln?PjUEe+|SYiyP-Ze{-(A* zYST#i8M&WeBXJ0;0w-9keS!89!^$0P9}KsEW9;QG1LM#(<6;oS9i*GIm)J-8v>nF@ z;Q>3o$Bw^>nCI?d*`MU@H@|P0-vM~X<3xIFuGmM62VvO$vxgXrJ2=g9GDq%rtVzo% zXgw^Wh`)MdNAZisvDLg)l@R)fRyO!J?+$BwV?!qql4u%jga*j0Kg^=FuCcb)h>+ZEy205BwMXwwFDck^+B2@ve9l-BXLU}%B?CM zO0F{tG_4DTXd;Hy6e?7NM#jj>Eah#6i|2$yB36f0ghP7>lr$RNZ-m3TgH)8JADs8n z-moW<#b|UW{^TSfBg|EdknqCE06VHP-6~Otp{Wx}-m#agl~5W=&*GF;nZblasSReE zkeyAAxmj~6ybAd&nx1eY%}sUjCHLk^+k%-_D?`df@J zW8|34n5`KG)UAw%B^YCNkQmF_bsm|DYwzAy-`BtRKjrGm5R@cgb1xQBZY`1Bi?gD zglB|^HbAfo%q;%o{$@ptsBIOHBwY+UF|uvVt3Bj_Pn!niM-NjR526at==kiRuc zirj;v^ruiqQDEOp2xqIE{7JM#f65%)hRf|BVnej>24|S_KU$aX#b3&5G_%^mXx&}@ zAaea6TK_Py3sFfGt9qA#5N0G|H`GSpzS7z-aCd1#82EWR%`qruP zbRT%z{FSVn$Ds@LyTj#|wODFBtTql{05Asjg9Xle#{4GW&VsCMCU&SjS^Xh{U9Y43 z3x;va!;xnV>)3FxQtcP)Z@-CM|BlF|U!mAntZX_TlTgejkKQ$oDSTt+;bYy8X9 z+csg!Ht83H@g}{_m+ROI2xu~CA5g~uetqpUg|>wG8|qX@tztU$!|v1+)kyS$g(eF` z-=@IBxkQgRM*b{e9%P5pN$zN^#x;cQekLaBG9uODn1^C2>A!JbqJF6?zl>8$v~^AjQ_5MPG`}!Zp`9y1~o^Wk6Ql z48_jCtIvj9B;^-$(e|POC6?%95+nJjL*or-d>5v_*76CW^Y~9fRef-d*-DBUXeV3A zEF0}(Avw%OjkCg#mJjbBKdpr+VxNX}qA-4ynsIpwZ9|B&gR@@Co=P0vQBI-ZyYVnJ zKPP`ee^|}+a9Ze;Z$WWvE^}YTWo}%P$!QEW^l}A#*Hx~YrYkF~IJZh>yQ807@Y%im%x2n^U3Bk=D=e3r zgEk$Kn@{~@y%cajEz<>wDCWB#jz3_|Tbzwu>@zUwpHhm&mZB057HVkyJ-CS1Vf>_O z|Lv8o|EyR3tgf|j$>@+~Fux2`*y+*=BDabo}?kcn^=a<71bwm{V9>8q4oM zyqoESo24qzZ6(}evxzE@WBAkgq620)$bu7lt)h#EaF^7>s=q-OiJR{i*73(L%)p{|vcSW2D>uVqDB{M0-SdxXdtq^L3}q4(>qt@FtXZ zQ$BaapuFD_8_UtFBfHL^hE2ANBGi7!XCg`%86U3?UvUT@)aasU1TN6~A-&sZt4L$} z1#4%e3U=aWyJ*hNLOtA_2|cQJPBNdcXn!-W`G=9? zG5D$8%Q?ATgma}WPH=I)Zho_9JOSEu88s+S5_W67c#$;z5$NJ9u^;tgaS7h(ZC}^* z&hI5&oW^E8a#WO!?!z4K{jvQd)`*$JtdI06n;~t6=f^a+wRHYU%ED@LdYCi=LQJw* zsqs4IGduB+@Y}iEu6u*+18}e@rB3@B zoGC&GhRQK`_}+=)V9AUXY^3ijZ{}MiIFCT7HSeH&T*z|G4bDrPAB=}~CVopux8%Q_ z_8N~TX6A#UID*}HY|x$FkQsy>R;=5D&j_dQ{m~e^JN-vR%NYPUuz`2dx=_B#PChOk zM1Arx99EQpjK=pwxFJEyv*IUoWBGJL(VWI`%p88`M$Td}UwFU$k1PS<5WArg_<((< z=kL<7D#xzG;7sN{yPOY6U<(RR?gP&{SD*CxD)WHbB#S;{nuOB^jjXP2ZiYFcNp*V? zP-nrTqSQyZU9z^o!`}ypm6L-fTKE>z&elpzt*D&xEdMxWd8Q}|;WIH?T2S_^<*Qjc zTh=+romTKPqhzP9*O`RIL-P~Ae%5kbgK&cD4r=llcY47cI3AwyL{)He{^1eS_!+YH z$;RhkvZ&3G_<1=JkA@b?hs5Ov{qWG5`2dNZH-f}B6dcZw)`ww!MB}3&9Eq=J2p=`% zy`|+@<0fFIwS1}>Li!~;q)DG9$)|iH))Y|Kozl|c+2ybrrke%Ypx9jEn z-@)-o+*G5U5U)ATO+Qrw;2Tb_vHUg(AuuK{H39KZQo~x8Zfz?W)z+Fm`CU7rF z^jxDP3Un%o0zFDHM4*=>diE)a0{u#&KpHlKQu0z7p*d#kZtL9SV5{o|v-0s%ZkNmQ zHSPLO!d13;#*77S2dxQ`~It$eTaI4`p#{;r7mEBoZe>h^jV%dJgP&Wp%%qE7=3{gu1 zLKZym5~O?BX)y5#F=CwbVUBN>3@^UjtQR)zt}L-rvxe*dS3kbpM;Nwe_w$6$*=e3( z>vx#eXWQU1KcMj~vsHHaCdg-&SSeg3HI|}SnYXlXT(H5$>b<<%Deq>Hvgycj?LT33 zc*Tar9ck{E_1C9)1ko~6SR62V1QQkJO3`}5^rr4L(_=qI8mtuf5QVi}X7~voMTFfs z-eKPV8(hw6WOVn3HPgDM*4zs4@fMP*>cx5#Hz5p|0&w=Keg)v;Rro%Cwgp454ctF} zyXBWg&e*~bXt~M2umZw-1`&pJj$v5YDskX)B)K4=vr94{n!2;)M=)x`ud)leI8#{F>~J|2CuhkCm*dZC5MWjC zE2!*`;#Yi+67wi1k~2pJgIzj;=bdmVvjB}?RQgnqhOaW>K374@uqsSbXOx}O_{{mfjZ3^19Q5&M~$Wj`~s>}Td(BpC;jnHjO4nPNYdBC_mf zW{*mN$;^z{&&(|QnVDrjGmm0Fw&6T@c{_shJPPO8>Xkpv!uaFae!-+amh<2O)Nr0! z!+EkS=Vn+=G?+6ui_`$YY#<0m&Pl#%g^qz)C?Wmt^dXpw9p8T}=8{>jb}p|A8OkKSwRuc_VpdswzivOa zo1Z|^$eQ+T{w8L7d%N?Rj{nUhO=o{w`e*a?G}gowVp8EYEoU()IkY5N#9i&&?BMRw zLe;s7F8cntPqh;&^t$?H)RcS&n$J|{SQgu*vB>h0+rI;SPZ zJk6)au4mQaCjUm_*)1L3@9|AFoGCz(OkV8$=@5c6{^WLZVZ`NIRpn1Raqd;k_%pBk zp!T|zUYrlvwv^hiY&j_MNNiPe?tL&1Tn&MujRwLRI4daFK%L7}9!Zcp^7#6LKIVs7 z_7PIqj76RHk;pxK(o(6`qF(^nQXoosma5_Gi*())$;vRw<6A>>he^r<^>4M1uhlY- zciN3#b6hJ>0aaUJM!0Ua?2Xy7^)WDx&winiF{x@Q&^B-qK8>0tDTA=e7e6&943<|( zRU7GT2{&6+RToPy@laxMtNjzmeoch==leIQ-bS1k;ms0ksjn>oAv(OyfYqxdmNg5Wx+$- z{nJo1vPl5Zh&+q}3mYcHx=-~@b1i86Mp;&JJqQB=?rm?7`|x-t>n49n3B90;7u~MQl}}v(ZH2LX_{o>_^iw!H#348JFrU z?xf;w9Nl*jUSar9V7?-cb5eOOck#$R7BNGDTJ!%SHnl(^cB6zPB&InKRr00QFVae6 zoOeX14AeT$=Q7SZMn$)Q=tVwNgY-OexvvQQLwGCjW!wLCiEz%AEE-zAJF@&?BbGt8 zcwcSfD)VxZfL%}QaSmaz&AJiufos&uBHV%?L+XqsYQB^yJrtg`YgMBwp|cxIJL6c& zZg`IimcpVi4z@ux>Iqa76Cv*PVBM^3Q|QK&9G0o$O$LC2wCH{j`=kFde1|;%<8=}o zDhz{CT6NbBVq@E9_gFAfd@FD z3g$A&xYba*P9jXrN*v6ra>Lt)4PG+c`Y~?TuI9UoYnUzx6%gWnpqAq!-4SPXhxJ*!c80Zid_j=f*zTE{Q{L~xTKPz3s`-I*2^E@i(a)lY zuDK1^O{)Jk+d8KUR}s)1v`2fwy7{7B%MzW5s5#-erD@kg=6r&L=S~|$y;DEUWVyk->O@INVAO0?5mO8ggNiQ#RQ@2&9M;bUWLYeY^PEI;7-)=e` zhdZ?E8pC{JKv{nc_r&csIh5&Uh2X038nX%p7nfCA zHN(}*7FVVlmS!>=&aKMJ-OalEGQQlE69r5FdcdBa%}V2>xiaPt-3^tAKQDRhUKHhW z1t=njC}7mP|3b5;lBq){L)=Pva3!hI|E2|=Hrj~PP{5~HKuWF#Uv6;x8|>-bVALx- zjQ@h3^xaDA_Y87#WE-Sn2gQ#7%&e@VL^+W~cj@G%5Ux2ehz~4lGf6M`nEX za83mHOb#HYTGa9R2vD2@)P{i42=HD=NA0$zy#i*tYlA>f$^Fh2*lJ_I}-0XXK1(B2pV9*zKK7cA@IP|^t$ScM zhyWcqz{(J?Gy+_o1N<-q+z|mT$^lk|fLkKKmvVrggaFNdjmN1u!0Hgt5CIBvfS-qe zizC3B?b)8+6#~wS06TMlyF)-t1o&+ZaBm2h5dnUl1Kbw^$|HcE1FQ)FAB_OF>q5XQ5#ZA~Ku-wR8388b01tK?L}R9N?)Ckct4`&jI>FKwAV@m;<~R0&a@{m*)WK5O8AzI41{q zEd*R00scG(7zzQGMu3Sqz#l@uMG@fLCE1=2hk#iT;Q1Wjy%10p0XF3TF>&Cacv=Lw zHwP#P0iTTkKh6OPL%^pZz@i*rd1h_f}kbfxZm>2=v9AIJy7#9Ji<^aZYB`D61 z03Xc(ri3~6e=8)Pyd0n;1Pn!hbXyiBr6J&z2(T>&C<_74M1cEqfHOkC;}PJeIY31S zcsK$q%>m8|0sj^OuFCCTZPx@wM_J<@u>m#{*!W(xt#8k6%dG5+tDe!ezcNifd?SR&GnsAuS&#gD^X16H zz}h@@A_VGt2s`j>edQJhZx8=a?Pvat0yyJ+o8qeU4c@Jpl|z(1-}>x^Z^o-K8>=!K zs+Js5cUzu9;5MI#uk)_lSHoA^Po26cQ~Z&?vTfWkj=J0 zyJEaaO#EimgH*wO&wMg#LHf<(P#1-%rON%y?I>cDEZt}dt!d*!R-@Ln#bOfG zmpv%OQu=~9Mic{8&hL5W-_Tc$e$04OU!tDKHPzOtOz(W_GaLRgzv*k~o30Y_Xrm5} z3x*V2a9(N#0j_}XW#;+TjooisFN16Vj5m6JLXBui&%xa19?YR`1PvqMY#=k;OWrj| z`nDaA4LN^RHn5s~2rri^b{+@ut-WTm(qpy{{Yq6w@^sYyRc6f zD8o39c}rAL)Y=F#BZtRwyV$WB6dL4EZm$9j1p4eimn>6Wno3k7Ie&XEDYG7h(Is%{ zyQ57}S^vc5$`EhbrZAU+)#jP=Q;cOoM_TDaADY7^XnbfE(lRo5+VljIcwy>TODhxz zYB;%|S2P1}Z8aqs{t;OgHrCXSSA;?~ujno06`{2(LTd#^h9B}_=UX9^tZRh8D@)JA zCVVBusAQD*H5EKd5b8T=Cf+vHceZI|>8lsvwHEi;Jw%{vtmK;->MsLrqit0GP(QM* zH=eX{q*w?N7wF2kKevSaAO{KEc-=Q{?)5p8`m;8K{uv|Z%E?WOa4`{i*ULH4?nZgL zLHZjx^9@~S^Dz;JkIpAIwY-9+gw0-yG_YV1g-{(k$fD5$e@$aFbJl>B9{He^*`-kg zc;q@{&kSJ|U&1IsYuOhHofl2~*0U;I2VDo(rl(}-~#1X#VYsdA?+(k9! zeh&Qrqqu{`3>Q~5cJe^$8C1|ah$AmNL_Q^jmD$Teb9jUf9cozL`zn|i_`QIEw<`UM z)do>~Y2n?bkb>I4Ng;oqf0eJY@STIk1m?{QTP2b8qkgMM)Nj?LdhG-btgb7){&HqU znH)YA6_mofVP0Er8$I8;U<0>%%QmF{1$dZ`;@GPmzwTBxs>M6`WgU|2MZxK7zG?iA z@&WHcZCr~V$blydZ0EoU170O?P<#QO8ZhrGc?Zv4Tc-HW!``^13o^y0X2W&d&kw_a z*>g{QIQ@$YML9nyz1PgAOO9o_d zuswX(y^61+SG~YF;M2)|>I`xt2R8dY}7x@dp?Z|(XzmU{PgJV)B4#ule zCst*;@E6iFnzkicZ~!K?CTd9N87;CgJ-MtKSbj_K(a{5?W*(7Ny-{*{6Luxa7vV)c zEs3leMp=K$Kuty0!F9%|=Y z0<#r^4}!s&Rnw<8ANRlV^iI14=P_8W0Bs-qgaQ==CM!@$AfcoJjZuvJ7+`JGeqe&N zd^2-?>61<4s{B~h^t@*1U|N-*-_<|1DwSX5=b1nrKoekA%Fc=z(tbW=Kp~t_=n(9V zpT%O|;qQO&_mBKNz+X3i+xV+l9gF>p=ac-s!r$Bcy~p1u;wSKT0)JEZJB`23@b^Xj z&g1W5{u=nZmcLv0yPdz-6$fwnPJo2qI0E<8es&(|q%-9$Of7i#H6 zn%j6yu4c_ z{_}{IEiN4VS^AAivDq}g*)cHq=TpYlm;MTcYOiS?D8N$5`(=9D^;*S(nm1WkZ(;$s zhavGENT-{G^O`qVe(yo7*$3nSmR~u&=+cpzFmpVs`!88N(!1(RCYG=fdQCKFULj@0~I-`FKiawd+aWL2?$=8r(`2)qL8zzvh)? zHJ{~?URF~<^D;E=tN6&@$LA@|b*G2s=+&)RP~^6L)KTK=<6{Lb+4&|XGk zJGLU}i%~c0^9{k=FM^jh)8Dz-tnMD)6`OKH*u62C;-V<{W~O*_6y&>B?_M7!_&xq_ zqu?(y#V>_HESrV#CF_KEmcjGBO5*I$nCB7SSF42EX{Z|T;hs=9^=l&+r+!<)8{gau z$8Xd=WW4YLh&_}`-=~3xGE7Qb|2a^J9QxA*$}Tr$oTzfy(Vc!}!nNG%up_ylL8fxk zFPg;r9az7$G8f-r;NN0+_uIeGJmB9ceBaCZ&U4H)d?$0{n!rNx9V&}CYlgxU9PsYY zC$lyWUKf3P>bB4E#h2W-`!K)bX!yAyOn=iS&Z0H?XrD3O3cr!XxF-fVsR4(-ybiF0 z|C^}~e?W&@akw@YF%H-3*jraIQs%o^I-g$1Ugln!1*+EQ5q<@$()#hn%4NUpB*x`v z>^mozaxa2U%q`)}@F4NzsA<%(2S9ysrnF2TwE- z)XRSXR2}#SmMJ+@lF+H01Yeh-#X80A*kqgY@|p#7oU!7d+nQ}~Ryo7_v=p0UfHZEKCinnt0w zZ@$Qwg#?Dk#Yy7pZ9n>+0T$pXvyY_qRY%p1JQUWx5`fuHsrsMqHr3mHP)CO$-~CP; zX|}`GTTn2Z-QuQ<5so^7jzS!igCW_^WfZ7X0M63HngOsvMxh|0ra}@IgXi&A1Lgl{ zwW4(-xCcZrd{@2FjBIj%Mpg){(Wd}ggKaJUmi7*Q4Tq;! z*Kx??ByoY8O(R#9H#AtwPY%;8u0sQC?v+J1r$c8{ee)OWF{3sj>!;%(<>;?vPHEih^EUC^~{;@Dk11w48 z3rkusf+Zcias*3iCIs`py-A|G&+T&N^lfvixYnH=+*w|A!y7Z~4!bkNz}%wEh1GKbm*<|5<)?!{tZv zqj`EC!H-@?%>UE)(chdYe)J=L4L`zO$xYcMc8afAr)<7&fh}A1>gkf#jP!|Z)f!ei z$JZynS^N0<_%ke;KQHYxOW_+>G@+}W@lC9nR>`VK?nIX@M1-`8+KeSMR&$bdiC(#Y zXr8O;!slvy@~i)F{4zHP!6qY&=)3rdc*_b3|X^H>I`EN~LIS9ZWLYiru5NjBlka;RhF-uIEPU zs|U+V+sIb+5OPM$EAIz>pK;ga`8|XmDqT&y(t2!8(M2wlFS2{7(0BnKkJbD$a{689 zVC^j7TMJ{Oy${d%FBjlqfoLg+r_uhnKnw7{|77+KK2PIEljjk>a!}jlTgyK|B?uQU zqfq9qb)|Z3n^J>XzPEtwC}`uv@V?T!NLUrj)4~G$Hc7_q8*(k@}ZIFD+JhMs$05_ zfMga8FjdlPP3|O%*XvD=bvEd2gWhzLd!ybq>TTFf?zV#EMgneSGU9-X zTSEGX^N%I0237ewr?Is8~wfW!}0dNm9$QvC;+Sq!518cK>aV|a`;aPnj7(JUA~EB4sgz#Tdex~q zm8!#u+XjHI0j%C-awH)NrxhEi_J@%d;E>;nK1$3yVu&MdEmGY83$x2P zU=e$;{en`%9WW=y6&Dm5vI$sLmOZ>H5*0mWNm0ZH&ABM3b#gj!dd@XQ0qx@(bF{Bh zt}$r8@l;XA*vfOpc!LPYECju!gu4VH;rLVa2Hry0fZ`@jrtgO_Xi#@R9d2G``ei1) znf@UkFbpB&4r*K@?uZIy!)=zY_+~e5?$d_6gxcy^`aqFF=+(!|d1L)!vivAyBwimcM(&fu1M8YDDpjCi;9$gib^Tsmt3Wjm9}&;_hEl7 z*cq1Qc9{~KVJX2`rz?T^t(&n+;TgN!c(@qJy@b{>Yx%bQl19Sin6VF6BLd+PMIc)#7^{7Q&G_V17tIP`z14+zr#M0@FK0uhs!n~E2{COCta zmpJx(J$k9XB=;}D$e)6_wK8-+TIzUHSUN2wnEOZFJ_u;x$;s=D`A*umXHGJ$W&}eX z;=~498Er5S4hX;ydZjVIuc={^7YkN&yncxx!(!macrP~?@j-}Ebd+g6lPotsR(vY2 zBLp_MF)SGFc-&ybk}3h7mly^avAg0kBU6{^R%tL*VZj6p!+ymtyi}_M+at|M)U#B4 zA%}%lqE`}qL$io~6>=Eld~vA(U4dT&9Szwzn&M5lLe!QKaJepAQ46&E6=~Ip zV)~sl=tOC>$v$bX$v)X#B*uy~jY!i7k02&##3Ub6@*Ab82*b$;IC^yR(QR+bA(XaI zCJ)3a!2)U#fTwU$~it{01*|WMbb2%=UCXEjh4AV^EPhQ`gchW z?ErhwMaJD(PP8C)Xp|`tmPS`HbN@W;;$xa7v`N#vCKx$inr3M(mpkc}>DW)Rq)iUB zZ`xflDoRB8hNX3i7AB9fwm@zLec=w2v7*Tiu?Y9!El{_si*>u27jGRRBK43*_y^qK zEt(Y>K~^IDlH8MKYl?dw-js&&M+N84grGe)jG;Fb4J*sT2sUIJo(sBR4Es(*v6RR? zzY*tbEt1fGLQPm35I5WvA{iCp^Hbp3GlIV%^^Do^euDc)Je*9cILSkMsrU;uhIS3< zdT%r7y0|3w3+Vb{G$TOakTujSWDQFgBN#*Q+B={r&>@Iq(;B?=9dWe>S#l; zCZrswl(TB|Fq%@`73-L~Oyvc|J9BuaI&*a8L)|{10E3pce4-Zv4b{)&dW*hp3d4qM zS(1}p#H-=s*)peD%FR~MZ|jPRlOI>9{KB9!^yM*xfcpM~#vnN~rH5J&9Y@aoO{h)I z7d@jzb!ZJed78}Ta3gYYbU9-1b8@!>8h}_=iO-ipUM?Xo8FR=vJgu8>dCyAN8*=p> z$OG(KL_kLr(=@b?RDzK!;$$H-Buz~&DdRrHKR`r5p=LB4!H>Tf#!x6=28}U0dWK7* zgZ$X-DMuy~`VgW?&NOSsQV(B<6@8+Sui~h#wkL(Lq2j$k2UNUANYFWkfbwM$q2WIU=69uzz7hv_>sH8x*4%1c;`nn(=f6WyyG58= zzh4I5@wJUX#R&S?r{z9MW3bLrxJ!J?_&i$u0Vi;S(R$I{`oruFU0}~+uFYUW?J6k$x8LG0#%DG_HYD=OejJ1P|HKHXfn+LN4q89c|+0Z3h5?q1xpGchR zCNOlBjIhA;0`vUtyB4Of-6R*v`#sL;Nw%~oRv{U|Kjjf4!?xKZ>|XGK4j5!`93jfC zpKF2)a3uMRc;g(hW`uHz8T;V*HOAiE*!d0=UK)li+1m}nJ$imY@1GS72L&wpv8hjv zce{gU^meOf^tJSQ4e}@Q^%3sf0LqX%7-0z6DOxW@{|s#F2+tcl@9?06HitUQ`GjKH z%Yz=EKL-omaz~#bWw|yNt=al>W#e)*UkgYGPzPE3?_s+nM#h5OsAfQ5MiVB>k&)jF zTQ3j-adDg%4UK)K21HQ^DjFMQ3@v7qklRdSA7~jY3p$2b2G(Nk$LurBaOOh4Hp3W@ zp$f}N952P~tYE`sHlr&y_XN7M`%uaw_dG0rSzaj7mV_7di{+0-%62iIz)aE$wAYWY zkoE#mUb_5KEF`|b2aSkMRp}QN_{DNiW(Y`rL76*sNuOBQLVO(}OzL7VlijH|d)%ol zUU%wkK5z^CC}Ha;EqKMf!5Jx^5J^?)sJgd?BOmDsxek}o z$kc$YO*yuFh%nsytQ%}lw&d&+0hcrAT0ks0{T6zmxXq9WB-zgzbsF2!Kyp8v?BS#8 zx(618_A|e{;k3D{|=t|N-u~mq|)#Rg~{zTBE)uNB71*>Cvmv zM&&4{aK3o8>k2J6y>6}CSQBq;)XHz93aF&#*i_8Q^FD2Kpm|N&7*~Gn@dk1cP(4Z9 z>qsMz(vxK4t%Z-*X6V&uT=EV0^L(kXJT1g1rE*vZ{-Z&vmk&OyuL{sbs;%B3%GY)u zD#r7*UU|_te*0k$1H0g1AMc{WUfLe#T+6w37a7b-*6fzhBrRLrMSY|`h2LzFGy1&G z9U|ud2OuCn4^xa*ayFz zI@s2+pNa#=i(GJZ&EVrC2TK>l0T7wYJPL`HQxfBzP;3&@cTh)^=5*>cZ+ zhSCNpilH5 zbU(IN^GMo&i2EgI!i%O!*Cre>BKK`Vko2d|nYTm>QR$g@;f@g1A589e>Z^%Lla$N_ z7roQ3Sz?@~xHFbk^c5sk^nwu0Fq$HUHw1}NXQ7G;Q5~)4nBUQ4VUUPSXKPAwd+ZqL z)|{xEFV}bCph8>+)UvpHW6f8PMx9)mMc*GM4Pz*_2MW>szov|+=!hIa9Z&{j32}@w z?wVz|B)LJ|UZF=5gXGvCc})+9CVz$V7zDsa8TBtZ`lRLq?xL))uX-?-;E^7kHiNyL!1D1JlmMv?B&{lE|dJq=?3 zEt7MnPF%k=osm|Ajd~OZdhg>}vzGP5!v zV^ffm@eT?N3}FtoX{kv18omg|C{27cbp{a`8Iit4=UCU=+t97iUr6_Tlp30ZR*}7^ zE4@dLvfNRTLiIT4n@3Q2quKpxP)6HK##9n9+9ugEl#QZl1zG5iL3V3NHSPHmT3Faf z5=j#_Bw5zdO7h+kHsI#6A9jI*&2lkuLv<(aTm<+FZmQj&OM@+K_z<)-^AFALV1WY* z&02GQmZ#_h*d4+bIgQq944x>2PGbau6fH$74OmZuFK2chvr{y)EO#DdXd{#39$8IT zE`y(5Y0u8%yq7!sOC4y3m#G4oEb7(8_A5E+<3D50l~`P~kZip|W-8JS^G;1WEIO;=HMC`gDTnH5 z?)d0Mrfn+%ur>xSjFp0-WVwr5IrIe1RwNQV!rpEejDd=bu)9ITRbR-sSw=QOm@?9l z30*z%VKPcln#AzK><%0-Fl1vOsE2kp6`}(4`hCOpv-JzbM2-Em-D2)~qDO;gdjagP z?a<&^;Q)b1!dBUcuD7+-fhdi<&CWKz?6hUSB@1XMrzv%@n0yAuV5HR>?0Q1JNQO~Y zejRx@_;*b6KD!+hCH)1&^)8qcMqb!o`<#E1-j-Tg(C_|kk&gQBTpT$oTRgb5 zyir#d=SSIM`_hT!XWe3J`HT$LuiqoXkBEw*g(EG{F@B@vVs z%J$l~*f%Yt5=kkg>Oz5tb4mr)r}WSo-HH)s(v=DCu$$ify=}yD!a-m*XHu6Ng*4bN ze#ayjb!h<_s)8?x8X-BXldzo-+=Iy=5a#fGl=X2`j|rG}J}Z1yVrHPlX7^WvUH9(9 z#%3#*UNPUtia9OnvKNXqoO=Yga;uO{#MtVJ5;2C`!LDbw5sX8H8Km#Nt0;j3Jxmnb z$k7NLK#S2C$09cNfuZ?Q6kg$voC*pS93ek1K&%`#I>kBB zI;|Z4m*kZlc+xne12*>Jl|)_l7 zD_3#y1;rI51N)gH>kB$UtekBA#C=yyYKpYDXz{>T<$x6j&E4`xpeGf;A6LrMh%_2lbx>==UyQT zlyU9l-f2!9%AUmBUM9`Z%C(14fn_`kdoR(ln@dFW{VUx!FxX2ik-=%I^x`fP+Y2^T zJbOgYb9t zHa*jMdF0@u;hA}26Ymi`F9YvELZ-JwapCvn%9V?m4K}l521t&_E@{%&cQZMG@3%`G ze+Zxf+`Vz09tfM1VkMK(0PAOT#|8m9-Y(5OQ!)UQhf{uJWhIN-jko8p(f-zB+kLHF+k9Vgfqoh~49%k0unkTRJr5gjW=7^0bDrDPx8R?S-N(A*r7|Ql` zIW#4w1axxY!Q}J?csi$hb!@WF^>!2qNx7hI%_*#fj5G%vnso@9ppmwhbjG}zeMB~P z$d}4AvcYXi3pnT=ubNVq!9LwvcIsbKHkVA5r;rEKH#`d z_@QvffWu&aa?+yU#@X&IUwk6DGuUfxkMz&FOZ!MOZ*EHQZjA{4mh4Om<$tKUFLq~2W5si(7P-vjHP zm2Nz~=l9@1`smqukKL&-f@Sm!E(7VOt}wRQiJTD zag2R{ICW^4yeEf7F!K-KNR z$Zu05-$XZeK0vwlKc)Z2EM`3bL# zl24G+R3hqCiEiq+n5RfDlCqMD`wB7d(^?VpenboqyAT3GkSX19xFGoApmetvk0RaK zDGsY}r%#p#3*<~fFS1=XwJ+kLmeaV_;|Ss%2d^W=5P=*dL`u3>8iSOr%yHi3-vcC#J^>MA;v3800?p=FgtMQ)0?wfZ;it>| zS~u+@o5ZI`ZkHJ{qF3@vr$Q5;UqG3maR>s<3|;GddO2oreqQ8RPU`E`tDyVpj}d2k ziA{e%SU109L`!`{sZ_2mY-5D^5sL4*|Kv|3!LV>0UG*bUL3%xj#VB{wfeqTa9A;J_ zT`XEU^?E(>$`wde%q6}Nmq!soJghV#EiFeT84}&zJb4LyNjct(KocS+hnqs8t)XwhI#y03kpSy$FQF+Y(<8zv=$WhX9 zB1LJjrd)O(y%aR`1B*IH>7RdkAN`3)X42iWt#n4fA;KfXFH)63gAj3!j_W}m!05Ow zqL~q?_&_3pzZ2{FA27;^I?w0BgbA3kqDTyM(7J2t()yt#dqZ-WFb;xu}3P zyf~f^R)vTd2sobY!x-m34lP{>S3(7hdrJetA7$*JN;||C;x`-ueW@8G1%GizZ#0Ta zde*qm#)i5NN+Q1xI0nSp3BD;x?!CHL! z`&xFZF!EtmTh(CTzaqOCmDg@oo<(-6GGsTYJRPmbZrNHHi$2P%g|;e-(Dp=!dg%MI_jKTgb`t{#1Qr_5Qp2^w6K>eW|=WNekqRedjf*U~o_Ws5OU>OFHPivra+b{lRFXD42EM^%weo(Ia0b8`(KxAsaCQ5+fTk{2D%! zJzJ;nvFyq4VwCVR zu#VEheJR;RY>G$6usgNXS2*xnPVm!FbEbpnF1Lcs#-k_<#p$(4b{kJ<&0rJJzNW&swRVxHn05Iaez zv<~i6dK7Qlzy_Ex*^4>Avm1K_ko}Hmq91T?kku{erERmG-whQgclmQcWw=G^v-wzxk}uEdzei)iX+Ed4qH7O3jI#%C*mnxD~~ak7U|cUxe}kA@$0i^ZAl4y@Y@UJC>jHER;ca8%9v1*pjYhwBxGiXVyws?XyunQ>_G{M~~ zQVE*^kCfKDCTb2~0ZgfA&F-a4lEisViEO--I-tMQYm`Q@GgGmR3?s6#7Yb z%z(v9#AcDl*%Ttzb{tK8fFiC`p*&D}7;{pKFRzt4M5~)IS&y^CgLrU82h?YWskshx z2}*whZ=!UHccP6>dKePa84@bMhwlGV^?e^WNWX|OLFvG8>oGDSeJpNA4YxXM&w(&&8D!CE?1QA7& zFtoC0`7hWN^3wS;0kv+nP_H-k_DWY>zjStTYJUBG={#kHP43vU{sNI}$bCXe4goeE3Al^* zj&SdMFDDR9*?1|vxw?N~01>Qnf6l&vDF>(0^G|t1zki|MkLvf+`u(bY|EqpKnJx1? zOX>N&=#4tz77ae%^C*SNpYdY81v6cA_)nq&>U#DQemW3NjJ z5od2Cp#rCcT~0Fohsio>@Ab#)$c;7r4LbhkaPKToZK<7O22SKZSu0t-lh=uNI|~NF ze;$lIlClN=g?b8x)$5)}!pUo}cR4-e>q?w-)dmKptR_O^dm{1Eg@b6rSr*MSq$mupAW}U_;4)|>B`V~D2_1f ze50`Qk7rn}DLkdvM{qbizA)ln+;C)20~`BT0#Z>Uu5jc_E&XBVdq}VBPoQIRl`_8;6)A|B7x#jfL&3}}9JVCR0%Agi&&Uy$P%lBU*85>=M2=mdkJi zC1vGezK8}IIT~xntbKH4#QN>J&9=r?;SYlNAs`lU@_LtMO27}s2-j*`qV2IMI77b4 z_vK;0QBD5_^I0aD#-{6XNIR;`^U?5m%y-SA!eOCoSP`@PpXT z2x2l@ZGTAuwdgp7{*FIMOyXDSLBs&9FVAOu?3?s98`wc`InbYX%D)NK*s;x6HXn_o zyi4(Rer$7ZxDdp~=i3$^NO&0(7QQ+C=UMvGHv(c4_2-6@b`Rio{B^@ajp1RAc=(9f z{{sp45}J}zCtUtt*>qZ)55{`H#N2#?@ZiqbD$KM>Pi&%5!sBRog-H4`AL|%M@A2aI zgq}M#{HBTD_ieh&5S7Nnx=ASZZC+mt_$fLH{O{5#ycOnxXvkV^v9jwnO{8ZwH$wWxfFDEUf0r=wF`4ge@*NQjoeQsjai}s{Lv{=5?h z`U}8w32S-Vc5z@t8tIczw2Cd_6Lgtma8cau((90auS|G!42{1%P?K^D1El$50ge~e zUs+a6nN;o`LzS$&vM{5qJea%J{V#=S?0;aNdM`%KDkOI>Dr0tw!A!y7a7JWtagZ@5D@?1!TJky&Rv6NriWf{ z8x<4Pts)vR!I@G8TDJ#DE8s3d?E$(?HA^#i!_zQ&_M@KLefsDC%S=$lq0TAGC1pkg z@B=+-;ul`H{=IGrm`iug!Eqq{M-z3mIwJ>Pyx7p;zMrUM^i$ejxVx8Z#{O048S6m$ zT0ApN6eA=vUq^97R-)w&q@N=$6vLcAXKnUuT%ReK6Xdh{TWVHFNBabVzwfK%Ep z2Cpbgze7j8u}2(f8)b}CBpJzwMi%bgCaKB>f-bEHFP0_^JH-Ee?s3$8cwDB8jg|a&A(K?PVlKlR`)A;xLwk7*%6IZ1J;W`0b2?t9IpMK0Q5%fU9gJ2}t0}20xNhM)83XCh>GgOS zTNZl;(z~%SAv06V&lu)tsn@VSE_N1UqO48s>GK9=VV=eIl39BKPtvsc5&MN$W(_3; zOJU+Pp<1#)zZVFc?!|hYSgH9N_i5pgOz&2HIEtpzHeuicQ>K#Ip+mcDw1;TIx4_ZF zO(-(UpGYpy>B^4%Ps6~WTH#Hzy01(~CU*Bd-Qs__8LmXbA)1e#_qB^iB|6s+I57^< zKuVth&`zPl&V4TCT3U`9LM!>KGPnb^@x&-_DHwK-!dFo=@M%H+?NcS&WTYta zreq|Z4fj)j9Nn5!JZ2AJ-=N(cNZ49rG_9hplqg2`Q-4ZxcmBOZ)Roc>;^{=pUFUtv za67kkEj}i-wi3Oqm*r}2U<@=;&Sq&Exp(rutM;{hGm@|n-<>_K0||K;gjsea_%Ico z+UE_=?(Fd%NXWzr|H0uLOuO;bTXf;kqJpDsFezfn?wixsmeQGWW+A-1>%8{?Qo%X2 zLxeFO^y8`AC-id?7uY8g51*t-={jNl15N*so(#o-_Pjx7y2F3 z?@P4*OZ9t#ejnBN$F%td{XV5JK5f2Azvt-pO8suv@5TCEtKWC&caDC)qr)2UOxKtU z{k~h{Z`bc@HKtpekLY)r#{3F19mL>oXnZu#H;A+tEbdcoBZeWTKxFe)%KvzIrLhQO z+BrX#aoi7X$ewfDXV|xi|LnN`=w`3jPmZwP6SqSYo@EA8dzbmDDlQzQLv$jj74EX> zoxB6-1LH(mruL1W!DWm*2BMoq@IC(cO@sOK(F%?~5VuD+&ldYVL-qh2!0}^4c4;-o z!$bCFiLZ9Z-XV5Ncpcx#8_o2hn|F?I_ltc)JbXmwkMi;y_kg3kW{bOX$UPwT)9($& zS0VOeL-E&%{gokmhuHTI*>!n8K4h12I{rZ19@#AM>-4vRBmG`+-w<~<>nA+Io1W67kA?q%`tV434!TD-yTopZZ7krpe**~)s z%`izi6?S@Q#qlgFt)C4`=yu<;C$a}Tx%uzk|28FX_dKY9-(zzB3cqb8ywijYCR}2| zc_zHXgkw#3hZ#?U2`fw(Frm+clTA3mgkQ`x(tXW@zc%4dOt{U2KR5k%oA4DA_M1>= zo>%Y+Bfj1%4Yk{K6(|G{fyO;TaRoHr=&nUcb7^2=^ls7MbuK6K*!) zArnqC;Y}vI)P(W)!XxY%at@N;A&{9;qk=%{b%9Mb8R~R9UG1KcAYTbvc_8)n{I9NhT3aF&Efi{ zMl)jXq_g{9QC;h;ZfIz#sWv=s|7-qQTB_H28$-8x8$yk%!qPvf@|I9&Ca}7;##dA4 zGvQ1V=9zGo3G+>OnF%j9;T0yFZNe)}I8$T{{WId9DM?xFc_ZwzM%eR5*e@GlzkG!K ziV=2^8T8Lg|H|R^yb<}!8*hi!{(^?+lZY>Y7TgyZ2 z*76X$wLHXbEf29<%R}te^1wc*sbLuc58i!-C7=E zx7L@~b$J9sE5p)GmejAR6T_umrSg0;FTZA1(TtX+rtplG>RV?B!gnR8nLc0MbU{lH z^Vi2AMEUg_u)D>tt*ZZD4LJ0U_l7cxm0#Q*H_N1xxD6P;Z^vP^YF71Xnp&+|dEbh; z5m=*|REr8>YE^5&)Tlbl4QeLuJRYB#u4=J2Nc8%53%*xI8wTisG$y`mwsaA`dJH>+<66*e?{$w%Bh=RAXNSDF1t9B+QdAm}SJ`TGI`>dqnV z`Ts`f@f-1vjIw{6i64J<_hA(?GVuQc_8*$%HZtlnWm*3@@V$-U>=6tciqCM+&7i~AXi zSprEM@t+WrnEgBwuk~lOZ6Zuyl1(8~SiKB0OS`MGxVK{-nYILwxZ<%`)0g^^* zlmX8rpOMbz1n_Y0-J9IZ_a0Mv6c9W`U%|H(xBKHCj~_{6dN|Y`&aSP;sC?kO{~pJ! z)*@vYVf=&sRySt*d@3G~5w@N5Y^tAlr`Ok3hpV%7-vFObvn8C*SF>`JQaR)G-GXZR z!X@R|6$=Y-Qoa#5LCyw8Hsb!4iNC;bf7c1FY#6RY$5$}|C&;Q1xP{>AhT{aKnup_b z8tubyL7hg&aGdtHX&A0d``bDkr~PdkhAY?pwu2k`-g5147r0&~k*Bvvc``dq8_gk) zhCRDn@+-VRZk{)BErNp=(vJ<|)-*MQSJz?}ekOd)^n2TFx7DtoN>oRdStl*wn(UH7 zou5s^aM$VlY_;Ob4S(B);S7K9HszG`!wszf^?=2mjUAq+T&8>Wbw-$7;I_F5({6lg zV5Ut$wx&HB;!-L<5Bd_EL0?){$)06XX+yYK;%_!M3mxVQPR5P@JbEy!Ay?8CD{z}` zB&m7AwzxgJMYmP&kX?U^rx3ov9Q^bZt=g{2SS^q`YBnROvZ*Qw)6&{flkL;;RTa1n zT}K)>C}ZL6Kl*jVJ!G;Dk(g)qwbQn{DXXgA@UjO~W|C7SKGRD0KL1Kcg|&v+dTnC&k^Sh+DvJlW=zuK2y>etato_?F=clUan`Wzv+pknp z>aJ2*wRfnj{B~7M+W<<} zT@-`!5;X@JHR#8wWIg*fo^qVB6b|_#aGK7X%7dRvM)%3}mxI|5^ z4XDZai)dO4R9gSIJ`Z#z<=rp^(rV*%Cyz5pj^zB9Tk8@%*f1wjHlk z5Aw24RJO@U#Cfqgw?DJ*oZj^A@y+9^Qp-HPtW1~6+~!fKn^IMB`#9yS8?UxLrqs(k zI|-VVk*G3yQj}-+SmoO3*7#3x>lU}m64hmAlGMeqF>3VDWF`3T{gQSnZscX8UE&-J zi@P-)+q@(-uV<9Hj79zPY*}hTUs~_DZu0l_aq8+flGK%(#wZ^f_;Tt}l=tL$ zDkFS8JmUi8D)%b?UUItZL8gu?o#RD(JCDtke)wErRSygS!t+;FjX(|vS+Bw zg`Q9zWsn>8sqESr>b(5vYC=t_Ek`BwyZRD)o!#~_Y4_?Po4UkCyC03=T%N=v4|P)L zoF{eCD6g^iI@B5a2T9JHM6YVMC#y_fO7Y~nQR@6->9n;pHKBIA8kav#ja@kUL%YdO zDu`+Y(dP*7t(JI;?JCP|p9;=)PLeu@_IXbJ7MCx1}Y)`^KuYC7%AVeaXE^-3HIJjZaYHC%gNEqm7|4W@Wfk26AEAD>mic4DnHsz|8_DfY|-86MUZLPYXyhdG+zd~JLyFq32 zn>E@@8aYWS=S;Tp#-^x@qf?cq=OX3WoufQADRnQ8$}#e9wP!9ETRf9`zvN_=y6~7+ zO$uK?U7jyxXPcyqI=P2%FT~64J}im%+0|UL?4*C=z;i;rO{Gp6qbB7~QIiNcyfloq zc9ftwZ}U1SaU=`t-Q!#!plXd&Y?94x&$2tC6B5@bKLSeXC2n;|U5=U>zFtkOy-rQd zuTWEMONPGLY{NbjhUZ+#tCJ@~_o+sA;551e=j?|ik-j`1lw`AOB5Ko|zwh7d+i`i@ z8FGk^x19@w0c4HVTsG&41hvj?uX8?Qf6$>`wcD>t`p4ww9b+eWegZ-NZd_j4TKj~F z89DX|DmX4Vnf#{F2h!K3(Z@<(DE#G+YpM86k{b1ENx-kGJb5?VV!CruVumWWr;ba2M@&-VY!g+YNk{Dxem$=T zc#@O4C49cb@?*!mU5F5O3#RuuA=Ojz8aG=_w>LU(N}OhQU7EBYd4ZbjvZWylLI0!s zM>Qu_CHS1go#nJE=dmfOH$s2LvzzF%lAS6UzUzeV3Qlc-Rs@fl4HxZE?9DuVg2Q>F zeOMTrt#MmsR^{4|yRSZcrKP*r2fgjU<}lCz%IPM%~J$99P?3Dub z)BSVQKZfV3X#PCqwar(lRh}}V58``j>3^VO>4Q>**M`$oMs1qP$RDpVK1}VC@2i-` z_f_*+&13ZS9lBddeBExJ;Dp1I_lXH=BK1cvKrc?+L0|dD&QX`XlA*l&Gs*WvHx{OC?dU=O~3Ag&#H~DfJ~kh^>||M%%uC{xCp)_%-^&MY=!qkLV9Q({1V< zY1^3?a(T?QJ<8G_b~)5Zd-MY5^u(#UKeSCqR1<77jF`HeWmz7V^1O1XN)AsW@6(kt z{|ePlTQAI2Y7S9Kxm)e27pABS>(bN(;hFR|KKh#J)aNC`Gws<;i>Z0&u<*m(9nfYEd zZea?woY?Q|Gy1sp8J0e7J(k4vc2(@u^1)Zp(@wl;dU{m+A0DKk%*EqESGJ$0^6Ji~ zeiu={#nf++uHShh>bGCgmHPb{!x&$@e)l9=>No08pV`+N_4^%ryV7zz?JP^`IgCJ5 zOIga{+Br()@=PKcsZ*=nW*LXt>~t^CwkO;mXGQ0vF;W=ipn|(~#XeZ4LYLIXe-$LY|7kQAL zHPrc`E!ol&Xv0sbblYytlWn7IYV->5E~j#}IGY_+cAqW9sZto{rpUNP=L>quPf+f0 zW15Lt$8UvV-bJ4-zdGhDM*bLLs!iZ`@f@>~3u%-)Ay<_ha9xOc(%Qf+Ft1@{@a4sbb^bWhG}?@`XZ!v4z8Rxj{(IkyQ4asR@0)pZ;{2;u z+#IS2vuN2=<1b$tZfdDs74rL+7M2wkR#bl7f1R43wl_k%Gb)lMD{Qhf- z3k&@L9Mx|U=Fgv1ZbTx%DkK~ZPp8lSS}0t;wmIZ4FJgXMn{c0B_|5#{Aw*&Ef@{Ng zs;qBZ+2r?2jG5ylePH0gUWiJJUbHd*wG~R(?H?Kn?XP7B;nnrq?z!s79yc zcTs% z)$hu|eK`KsaIN3Jys+{c71yahIOn%Ew}or`^Rrh}x7OEGwuW0m)vH&ut@IZyo1a}% zSKU$>ZmF&hxB82WfW}*L7bueS-Zo-r{IA4NxmZ51a+x+>L$II_; z4p~@Sc&%@yzwGit(_{kwby6D6|0{87I;eKer(*x}UnJR>nMZAti#o^rl?`pJb^qm9 zD4oac^J|*Oh^D$(*Z!^LM1$cN@X~7+md(uf8B~9vs)RgiE5o6d)%A_g39_Dn8KmvX znx^KpEF*jgUmK!Eo7N&FRW59AFRyE9y0zH4F|AV2i(0BzhXU1&wGE+`d=Vm7hgR1} zeDlNvbt2*u^TNv3=1@)jN?J>u;aRD3!Pm64w1gThElo{283CYq3f)wz7UzTT7_2vATs4 zs;q8VwT4qo3Kqw%k?Mx}+XkINjce*#ni^O0!N~BcauxBLzT%;+UM@G{kg8Kx44ErK zY8hWsPaM>)Hk2EAUAelcmale_)8@)b$w`eI1R1oEuLE7 z9jvUjQnnN}NHJ7GG9*#y1;Ywa<{mE2TwQR_aVZGm?jzrqxEC@`sg2 zQ=^VFEW2rhW`@!?T~@+^^_Ir-$z=<*pmU*CsUbS9ZmyrUx^>#E^^Mc;F|DPoFzuYbyg$XrB?r%JI$!R$?9Q5c%vk&9#ofF*wEB!HV-x1Rtyoh zENp=;D9@$z!+MM7en&7=y(aV}{4I`>K@u9_SlC#*w7I@iVRdtuR=`#$-A;cexL{M$Ep5$<%(~>`Ti^?S z6qjr3Tf%MC4a*v-zFN*UI!aoLLMz%g2FYAg%Obvqe8pq+O70gLOvm)2ew0ui=ZE!48MEYz~HiT2)D!b@E7d$*-_eWd0r@V0Q<=BS zz-(QgkadeTo;ubb1fJlz9kX>#Q|3C~yp(TYo>AbtJPnxpffHv?w#4ZLw$0?*`Nj_5 ztFsJqAMobOh#!12u;mKk!G8zPb`^U%NXrAPD`0M%{9EVUWIoO_mvUfEPT&zuO8M=xE% zJQR3=Nfk={9J4@{&QvdEKkxw011!t-0y~ydr{qDP%vH5d#teLxCmFNMb4_DzO3X4x zbvKXf1=pD)#+`2YnUAM8W>Bw5os3`4$^&Wv2@FY{70&93CVHS9V=TgjlV!n?h z2g=#He@FK9JjzoJUf^+_I?VmR$A1iez}yX-e?R4p*$=#O2ThuE{lJ$WU=0v+AF%!* zzL7d=27dMk-$FU4f2PhKGtAa~2G%_Y*8K^xZ(-$Agog$MUh-?&FL4TdsGD|&S>T+# z>xrV2DTj)Bz8g zc-e<_+{C{J{Mw(1i#TPU!vUUJ%zeP={j@#Ie&9_!n=s2>&q3b!YjsIrH0B zH+a_=bAod(nmOHfkvG&>(^SjwVeUoC$`?(WeUX=GsK(mr21ammFIpRFy=eZtl%zS# z>V;OXXjtndg2vXl7qzuCUe#Jt7g}B2I&F1*O-oa2)5`EPhKW~Ix2~SPX68j+gzfs3 z%siMgkVa}m<@L@9x3sl}OPDz{W1IXR#5QZP4oAq=5EBjc;kBk6OiSonZKM*aEo-S? zgH*mM)H>*BbuMlvsEpHsp*5ifZ-cz&UR2#$(zvGSmQc$@-nRO}8X0QLy=Y~1Lu=?F z?~K8C=FIr&ROZYWmZv#020xDAb7sWL4a@LKJiE{oL+e6zQ>VBx_p`3^P4+kFZeK_`T|HEe< sRzGw7tm9{!ezx^z+kPe|64QvYV*94;yS5+Q-oO3qN&h$Y|NRp9Uyp!|UH||9 literal 0 HcmV?d00001 diff --git a/include/GaussNewton.h b/include/GaussNewton.h new file mode 100644 index 0000000..ce17b48 --- /dev/null +++ b/include/GaussNewton.h @@ -0,0 +1,89 @@ +#pragma once + +#ifndef GNM_h +#define GNM_h + +#include "utils.h" + +#include + +class GaussNewton { +public: + GaussNewton(int L_, double* parma_, double (*Func_)(double, double*), + double* (*Gunc_)(double, double*)); + ~GaussNewton(); + +public: + void addData(double x, double y); + double* solve(); + +public: + double eps = 1e-5; + double* parma; + double (*Func)(double, double*); + double* (*Gunc)(double, double*); + int L, maxIter = 10; + std::vector data; + +private: + void calmJ_vF(); + void calmH_vG(); + +private: + Eigen::MatrixXd mJ; // 雅克比矩阵 + Eigen::MatrixXd mH; // H矩阵 + Eigen::VectorXd vF; // 误差向量 + Eigen::Vector3d vG; // 反馈向量 +}; + +GaussNewton::GaussNewton(int L_, double* parma_, double (*Func_)(double, double*), + double* (*Gunc_)(double, double*)) { + L = L_; + parma = parma_; + Func = Func_; + Gunc = Gunc_; +} + +GaussNewton::~GaussNewton() {} + +void GaussNewton::addData(double x, double y) { data.push_back(Eigen::Vector2d(x, y)); } + +void GaussNewton::calmJ_vF() { + double x, y; + double* resJ; + + mJ.resize(data.size(), L); + vF.resize(data.size()); + + for (int i = 0; i < data.size(); i++) { + Eigen::Vector2d& point = data.at(i); + x = point(0), y = point(1); + resJ = (*Gunc)(x, parma); + for (int j = 0; j < L; j++) mJ(i, j) = resJ[j]; + vF(i) = y - (*Func)(x, parma); + } +} + +void GaussNewton::calmH_vG() { + mH = mJ.transpose() * mJ; + vG = -mJ.transpose() * vF; +} + +double* GaussNewton::solve() { + Eigen::VectorXd vX; + vX.resize(L); + + for (int k = 0; k < maxIter; k++) { + calmJ_vF(); + calmH_vG(); + + vX = mH.ldlt().solve(vG); + if (vX.norm() <= eps) return parma; + + for (int i = 0; i < L; i++) parma[i] += vX[i]; + } + + return parma; +} + +#endif \ No newline at end of file diff --git a/include/LevenbergMarquardt.h b/include/LevenbergMarquardt.h new file mode 100644 index 0000000..25bbbb8 --- /dev/null +++ b/include/LevenbergMarquardt.h @@ -0,0 +1,121 @@ +#pragma once + +#ifndef LMM_h +#define LMM_h + +#include "utils.h" + +#include +#include +#include + + +class LevenbergMarquardt { +public: + LevenbergMarquardt(int L_, double* parma_, double (*Func_)(double, double*), + double* (*Gunc_)(double, double*), int type_ = 0); + ~LevenbergMarquardt(); + +public: + void addData(double x, double y); + double* solve(); + +public: + double mu = 1., eps = 1e-5; + double* parma; + double (*Func)(double, double*); + double* (*Gunc)(double, double*); + int L, type_, maxIter = 20; + std::vector data; + +private: + void calmJ_vF(); + void calmH_vG(); + double calMse(double* parma_); + +private: + Eigen::MatrixXd mJ; // 雅克比矩阵 + Eigen::MatrixXd mH; // H矩阵 + Eigen::VectorXd vF; // 误差向量 + Eigen::VectorXd vG; // 左侧 +}; + +LevenbergMarquardt::LevenbergMarquardt(int L_, double* parma_, double (*Func_)(double, double*), + double* (*Gunc_)(double, double*), int type_) { + L = L_; + parma = parma_; + Func = Func_; + Gunc = Gunc_; + this->type_ = type_; +} + +LevenbergMarquardt::~LevenbergMarquardt() {} + +void LevenbergMarquardt::addData(double x, double y) { data.push_back(Eigen::Vector2d(x, y)); } + +void LevenbergMarquardt::calmJ_vF() { + double x, y; + double* resJ; + + mJ.resize(data.size(), L); + vF.resize(data.size()); + + for (int i = 0; i < data.size(); i++) { + Eigen::Vector2d& point = data.at(i); + x = point(0), y = point(1); + resJ = (*Gunc)(x, parma); + for (int j = 0; j < L; j++) mJ(i, j) = resJ[j]; + vF(i) = y - (*Func)(x, parma); + } +} + +void LevenbergMarquardt::calmH_vG() { + mH = mJ.transpose() * mJ; + vG = -mJ.transpose() * vF; +} + +double LevenbergMarquardt::calMse(double* parma_) { + int n = data.size(); + double x, y, mse = 0; + + for (int i = 0; i < n; i++) { + Eigen::Vector2d& point = data.at(i); + x = point(0), y = point(1); + mse += pow(y - (*Func)(x, parma_), 2); + } + + return mse / n; +} + +double* LevenbergMarquardt::solve() { + double v = 2, cost; + double* parmaNew; + + Eigen::VectorXd vX; + Eigen::MatrixXd mT, mD = Eigen::MatrixXd::Zero(); + vX.resize(L); + mD.resize(L, L); + for (int i = 0; i < L; i++) mD(i, i) = 1; + + for (int k = 0; k < maxIter; k++) { + calmJ_vF(); + calmH_vG(); + + if (type_ == 1) { + mT = mJ * mJ.transpose(); + for (int i = 0; i < L; i++) mD(i, i) = sqrt(mT(i, i)); + } + + mH = mH + mu * mD; + vX = mH.ldlt().solve(vG); + if (vX.norm() <= eps) return parma; + + for (int i = 0; i < L; i++) parmaNew[i] = parma[i] + vX[i]; + cost = calMse(parma) - calMse(parmaNew); + cost /= vX.transpose() * (mu * vX + vG); + } + + return parma; +} + +#endif \ No newline at end of file diff --git a/include/clip.h b/include/clip.h new file mode 100644 index 0000000..96bf93c --- /dev/null +++ b/include/clip.h @@ -0,0 +1,92 @@ +#pragma once + +#ifndef clip_h +#define clip_h + +#include + +#define INF 1e9 + +class SigmaClip { +private: + int maxiters; + double sigma, minValue, maxValue; + double (*cenF)(double *, int) = nullptr; + double (*stdF)(double *, int) = nullptr; + +public: + SigmaClip(double sigma = 3, int maxiters = 5, double (*cenF)(double *, int) = nullptr, + double (*stdF)(double *, int) = nullptr); + ~SigmaClip(); + double *clip(double *data, int L); + +private: + void computeBound(double *data, int L); +}; + +double GetAvg(double *data, int L) { + double sum = 0; + for (int i = 0; i < L; i++) sum += data[i]; + return sum / L; +} + +double GetStd(double *data, int L) { + double sum = 0; + double mean = GetAvg(data, L); + for (int i = 0; i < L; i++) sum += (data[i] - mean) * (data[i] - mean); + return sqrt(sum / L); +} + +SigmaClip::SigmaClip(double sigma, int maxiters, double (*cenF)(double *, int), double (*stdF)(double *, int)) { + this->sigma = sigma; + this->maxiters = maxiters; + this->minValue = INF; + this->maxValue = -INF; + + this->cenF = cenF == nullptr ? GetAvg : cenF; + this->stdF = stdF == nullptr ? GetStd : stdF; +} + +SigmaClip::~SigmaClip() {} + +void SigmaClip::computeBound(double *data, int L) { + double std = (*stdF)(data, L); + double mean = (*cenF)(data, L); + + minValue = mean - (std * sigma); + maxValue = mean + (std * sigma); +} + +double *SigmaClip::clip(double *data, int L) { + int tmp, cnt; + double *value; + + for (int i = 1; i <= maxiters; i++) { + computeBound(data, L); + + cnt = 0; + for (int j = 0; j < L; j++) + if (data[j] >= minValue && data[j] <= maxValue) cnt += 1; + + tmp = L, L = cnt, cnt = 0; + value = new double[L]; + for (int j = 0; j < tmp; j++) + if (data[j] >= minValue && data[j] <= maxValue) { + value[cnt] = data[j]; + cnt += 1; + } + + delete data; + data = new double[L]; + for (int j = 0; j < L; j++) data[j] = value[j]; + delete[] value; + } + + value = new double[L + 1]; + value[0] = L; + for (int j = 0; j < L; j++) value[j + 1] = data[j]; + + return value; +} + +#endif diff --git a/include/getADC.h b/include/getADC.h new file mode 100644 index 0000000..80d778a --- /dev/null +++ b/include/getADC.h @@ -0,0 +1,192 @@ +#include "LevenbergMarquardt.h" +#include "utils.h" +#include +#include +#include +#include +#include + +#include +#include + +using namespace std; +using Eigen::MatrixXd; +using Eigen::VectorXd; + +void getADC(const char *fin) { + // read file + TFile *fRun = new TFile(fin); + TTree *t = (TTree *)fRun->Get("Tree1"); + + // data numbers + Long64_t ntot = t->GetEntriesFast(); + + // five X-4 block + UInt_t bArray[16]; + UInt_t cArray[16]; + UInt_t dArray[16]; + UInt_t eArray[16]; + UInt_t fArray[16]; + + // calibration coefficient + Float_t CalC[5][8]; + Float_t Calk1[5][8]; + Float_t Calk2[5][8]; + + // read adc data + t->SetBranchAddress("adc0ch0", &bArray[0]); + t->SetBranchAddress("adc0ch1", &bArray[1]); + t->SetBranchAddress("adc0ch2", &bArray[2]); + t->SetBranchAddress("adc0ch3", &bArray[3]); + t->SetBranchAddress("adc0ch4", &bArray[4]); + t->SetBranchAddress("adc0ch5", &bArray[5]); + t->SetBranchAddress("adc0ch6", &bArray[6]); + t->SetBranchAddress("adc0ch7", &bArray[7]); + t->SetBranchAddress("adc0ch8", &bArray[8]); + t->SetBranchAddress("adc0ch9", &bArray[9]); + t->SetBranchAddress("adc0ch10", &bArray[10]); + t->SetBranchAddress("adc0ch11", &bArray[11]); + t->SetBranchAddress("adc0ch12", &bArray[12]); + t->SetBranchAddress("adc0ch13", &bArray[13]); + t->SetBranchAddress("adc0ch14", &bArray[14]); + t->SetBranchAddress("adc0ch15", &bArray[15]); + + t->SetBranchAddress("adc0ch16", &cArray[0]); + t->SetBranchAddress("adc0ch17", &cArray[1]); + t->SetBranchAddress("adc0ch18", &cArray[2]); + t->SetBranchAddress("adc0ch19", &cArray[3]); + t->SetBranchAddress("adc0ch20", &cArray[4]); + t->SetBranchAddress("adc0ch21", &cArray[5]); + t->SetBranchAddress("adc0ch22", &cArray[6]); + t->SetBranchAddress("adc0ch23", &cArray[7]); + t->SetBranchAddress("adc0ch24", &cArray[8]); + t->SetBranchAddress("adc0ch25", &cArray[9]); + t->SetBranchAddress("adc0ch26", &cArray[10]); + t->SetBranchAddress("adc0ch27", &cArray[11]); + t->SetBranchAddress("adc0ch28", &cArray[12]); + t->SetBranchAddress("adc0ch29", &cArray[13]); + t->SetBranchAddress("adc0ch30", &cArray[14]); + t->SetBranchAddress("adc0ch31", &cArray[15]); + + t->SetBranchAddress("adc1ch0", &dArray[0]); + t->SetBranchAddress("adc1ch1", &dArray[1]); + t->SetBranchAddress("adc1ch2", &dArray[2]); + t->SetBranchAddress("adc1ch3", &dArray[3]); + t->SetBranchAddress("adc1ch4", &dArray[4]); + t->SetBranchAddress("adc1ch5", &dArray[5]); + t->SetBranchAddress("adc1ch6", &dArray[6]); + t->SetBranchAddress("adc1ch7", &dArray[7]); + t->SetBranchAddress("adc1ch8", &dArray[8]); + t->SetBranchAddress("adc1ch9", &dArray[9]); + t->SetBranchAddress("adc1ch10", &dArray[10]); + t->SetBranchAddress("adc1ch11", &dArray[11]); + t->SetBranchAddress("adc1ch12", &dArray[12]); + t->SetBranchAddress("adc1ch13", &dArray[13]); + t->SetBranchAddress("adc1ch14", &dArray[14]); + t->SetBranchAddress("adc1ch15", &dArray[15]); + + t->SetBranchAddress("adc1ch16", &eArray[0]); + t->SetBranchAddress("adc1ch17", &eArray[1]); + t->SetBranchAddress("adc1ch18", &eArray[2]); + t->SetBranchAddress("adc1ch19", &eArray[3]); + t->SetBranchAddress("adc1ch20", &eArray[4]); + t->SetBranchAddress("adc1ch21", &eArray[5]); + t->SetBranchAddress("adc1ch22", &eArray[6]); + t->SetBranchAddress("adc1ch23", &eArray[7]); + t->SetBranchAddress("adc1ch24", &eArray[8]); + t->SetBranchAddress("adc1ch25", &eArray[9]); + t->SetBranchAddress("adc1ch26", &eArray[10]); + t->SetBranchAddress("adc1ch27", &eArray[11]); + t->SetBranchAddress("adc1ch28", &eArray[12]); + t->SetBranchAddress("adc1ch29", &eArray[13]); + t->SetBranchAddress("adc1ch30", &eArray[14]); + t->SetBranchAddress("adc1ch31", &eArray[15]); + + t->SetBranchAddress("adc2ch0", &fArray[0]); + t->SetBranchAddress("adc2ch1", &fArray[1]); + t->SetBranchAddress("adc2ch2", &fArray[2]); + t->SetBranchAddress("adc2ch3", &fArray[3]); + t->SetBranchAddress("adc2ch4", &fArray[4]); + t->SetBranchAddress("adc2ch5", &fArray[5]); + t->SetBranchAddress("adc2ch6", &fArray[6]); + t->SetBranchAddress("adc2ch7", &fArray[7]); + t->SetBranchAddress("adc2ch8", &fArray[8]); + t->SetBranchAddress("adc2ch9", &fArray[9]); + t->SetBranchAddress("adc2ch10", &fArray[10]); + t->SetBranchAddress("adc2ch11", &fArray[11]); + t->SetBranchAddress("adc2ch12", &fArray[12]); + t->SetBranchAddress("adc2ch13", &fArray[13]); + t->SetBranchAddress("adc2ch14", &fArray[14]); + t->SetBranchAddress("adc2ch15", &fArray[15]); + + TH1F *Left = new TH1F("Left", "Left", 300, 0, 300); + TH1F *Right = new TH1F("Right", "Right", 300, 0, 300); + + for (int i = 0; i < ntot; i++) { + t->GetEntry(i); + Left->Fill(bArray[0]); + Right->Fill(bArray[8]); + } + + // use matrix and least square method to get coefficient + // SigmaClip *clip = new SigmaClip(5); + int cnt = 0; + double vX[300]; + for (int i = 2; i < 300; i++) + if (Left->GetBinContent(i) > 0) { + vX[cnt] = i; + cnt += 1; + } + + VectorXd mY(cnt); + MatrixXd mX(cnt, 3); + MatrixXd mB(3, 1); + double n; + double *parma = new double[3]; + double sigma, b0, b1, b2; + + cnt = 0; + for (int i = 2; i < 300; i++) { + n = Left->GetBinContent(i); + if (n == 0) continue; + mY(cnt, 0) = log(n); + mX(cnt, 0) = 1, mX(cnt, 1) = i, mX(cnt, 2) = i * i; + cnt += 1; + } + + mB = (mX.transpose() * mX).inverse() * mX.transpose() * mY; + b0 = mB(0, 0), b1 = mB(1, 0), b2 = mB(2, 0); + sigma = -1 / b2; + parma[1] = b1 * sigma / 2; + parma[0] = exp(b0 + parma[1] * parma[1] / sigma); + parma[2] = sqrt(sigma); + + cout << parma[0] << ", " << parma[1] << ", " << parma[2] << endl; + + LevenbergMarquardt LM(3, parma, Gaussian, GaussianJacobian); + for (int i = 2; i < 300; i++) { + n = Left->GetBinContent(i); + if (n == 0) continue; + LM.addData(i, n); + } + parma = LM.solve(); + cout << parma[0] << ", " << parma[1] << ", " << parma[2] << endl; + + // TCanvas *c1 = new TCanvas("c1", "Fitting with Gaussian function"); + // c1->Divide(1, 2); + + // c1->cd(1); + // c1->SetGrid(); + // Left->Draw(); + + // c1->cd(2); + // c1->SetGrid(); + // Right->Draw(); + + // TF1 *GausFunc = new TF1("GausFunc", "[0]*exp(-(x-[1])*(x-[1])/2.0/[2]/[2])", 0, 300); + // GausFunc->SetParLimits(0, 1000, 2400); + // GausFunc->SetParLimits(1, 100, 200); + // GausFunc->SetParLimits(2, 0, 15); + // Left->Fit(GausFunc, "R"); + // Right->Fit(GausFunc, "R"); +} diff --git a/include/utils.h b/include/utils.h new file mode 100644 index 0000000..06f3942 --- /dev/null +++ b/include/utils.h @@ -0,0 +1,20 @@ +#pragma once + +#ifndef utils_h +#define utils_h + +#include + +double Gaussian(double x, double* p) { + return p[0] * exp(-(x - p[1]) * (x - p[1]) / (2 * p[2] * p[2])); +} + +double* GaussianJacobian(double x, double *p) { + double* resJ = new double[3]; + resJ[0] = -Gaussian(x, p) / p[0]; + resJ[1] = -(x - p[1]) * Gaussian(x, p) / (p[2] * p[2]); + resJ[2] = -(x - p[1]) * (x - p[1]) * Gaussian(x, p) / (p[2] * p[2] * p[2]); + return resJ; +} + +#endif diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..02eeceb --- /dev/null +++ b/main.cpp @@ -0,0 +1,6 @@ +#include "getADC.h" + +int main() { + getADC("2016Q3D/root/raw/201609Q3D1002.root"); + return 0; +} diff --git a/process.cpp b/process.cpp new file mode 100644 index 0000000..5df244d --- /dev/null +++ b/process.cpp @@ -0,0 +1,9 @@ +#include +#include +#include + +void process() { + gInterpreter->AddIncludePath("./include"); + gROOT->ProcessLine(".L getADC.cpp"); + gROOT->ProcessLine("getADC(\"2016Q3D/root/raw/201609Q3D1002.root\")"); +} \ No newline at end of file diff --git a/readme.md b/readme.md new file mode 100644 index 0000000..e69de29