From 85abe25401c5cc330fec3fc500c9aa2f365876a0 Mon Sep 17 00:00:00 2001 From: Cameron Clarke Date: Sun, 28 Jun 2020 18:23:42 -0400 Subject: [PATCH] Updating Qsim with better messenger commands, new macros, and a wrapper script and analysis script (all in macros). Including the qe.txt and UVS text files allows for clean reproduction of the results in https://arxiv.org/pdf/1710.07100.pdf Detector construction needed to be reworked a bit to get that to work, but I can reproduce the 2015/6 results nicely now. --- UVS_45total.txt | 2046 +++++++++++++++++++++++++++ include/qsimDetectorConstruction.hh | 40 +- include/qsimMessenger.hh | 6 + include/qsimOpticalPhysics.hh | 10 +- macros/get_pe.C | 351 +++++ macros/preserve.mac | 31 + macros/qsim.sh | 114 ++ macros/vis_brad.mac | 23 + qe.txt | 548 +++++++ qsim.cc | 9 +- src/qsimDetectorConstruction.cc | 1971 ++++++++++++++++++++------ src/qsimMessenger.cc | 26 + src/qsimOpticalPhysics.cc | 29 +- src/qsimPrimaryGeneratorAction.cc | 33 +- 14 files changed, 4789 insertions(+), 448 deletions(-) create mode 100644 UVS_45total.txt create mode 100644 macros/get_pe.C create mode 100644 macros/preserve.mac create mode 100755 macros/qsim.sh create mode 100644 macros/vis_brad.mac create mode 100644 qe.txt diff --git a/UVS_45total.txt b/UVS_45total.txt new file mode 100644 index 0000000..90c7015 --- /dev/null +++ b/UVS_45total.txt @@ -0,0 +1,2046 @@ +178.41 0.618527 0.427469 +178.791 0.622443 0.427689 +179.172 0.626332 0.427909 +179.553 0.174379 0.562368 +179.934 0.628751 0.654833 +180.316 0.473538 0.489524 +180.697 0.559314 0.670673 +181.078 0.734451 0.743916 +181.459 0.698905 0.728918 +181.84 0.547226 0.592743 +182.602 0.222801 8.88282 +183.364 0.301337 0.458317 +183.744 0.532319 0.236116 +184.125 0.656703 0.23355 +184.506 0.616146 0.213931 +184.887 0.58486 0.0833822 +185.268 0.553837 0.132498 +185.648 0.593833 0.0932494 +186.029 0.57657 0.0970726 +186.41 0.66038 0.069334 +186.791 0.714758 0.0818225 +187.171 0.567066 0.0868458 +187.552 0.475665 0.0475853 +187.932 0.530196 0.0511289 +188.313 0.47674 0.0355922 +188.694 0.529922 0.0449169 +189.074 0.529284 0.0394132 +189.455 0.509442 0.0467865 +189.835 0.561462 0.0576489 +190.216 0.551772 0.041978 +190.596 0.517079 0.0389335 +190.977 0.510455 0.0475364 +191.357 0.53227 0.0362066 +191.737 0.507464 0.0339295 +192.118 0.5049 0.0387569 +192.498 0.484586 0.0452928 +192.878 0.499961 0.0328762 +193.258 0.48113 0.0335352 +193.639 0.487481 0.0362805 +194.019 0.491724 0.037861 +194.399 0.501482 0.04539 +194.779 0.472907 0.0339923 +195.159 0.506097 0.045871 +195.54 0.472705 0.0347932 +195.92 0.473513 0.034049 +196.3 0.463116 0.0378648 +196.68 0.487493 0.0526724 +197.06 0.474403 0.0469343 +197.44 0.471436 0.0437746 +197.82 0.465097 0.0364796 +198.2 0.449774 0.0376481 +198.58 0.43222 0.0350726 +198.959 0.4537 0.035721 +199.339 0.45174 0.037052 +199.719 0.439143 0.031043 +200.099 0.449441 0.0430499 +200.479 0.440087 0.0349413 +200.858 0.418757 0.0343527 +201.238 0.41468 0.0372159 +201.618 0.441878 0.0338896 +201.998 0.40779 0.0435623 +202.377 0.428405 0.0373526 +202.757 0.411991 0.0418579 +203.136 0.407499 0.0475428 +203.516 0.413722 0.034678 +203.896 0.394132 0.0366363 +204.275 0.383276 0.039823 +204.655 0.394775 0.035215 +205.034 0.374998 0.0432391 +205.414 0.380584 0.0317528 +205.793 0.370382 0.0422066 +206.172 0.375048 0.0428886 +206.552 0.358232 0.0287338 +206.931 0.383761 0.0343891 +207.31 0.367313 0.0348063 +207.69 0.364296 0.0312534 +208.069 0.364716 0.0330064 +208.448 0.361802 0.0374477 +208.827 0.373072 0.0327706 +209.207 0.374957 0.0370455 +209.586 0.361841 0.0392752 +209.965 0.356173 0.0387103 +210.344 0.351553 0.0357664 +210.723 0.340085 0.0316846 +211.102 0.334714 0.0309992 +211.481 0.322511 0.0311813 +211.86 0.330135 0.0319701 +212.239 0.335948 0.034297 +212.618 0.330969 0.0348467 +212.997 0.321273 0.0400792 +213.376 0.318524 0.0396558 +213.755 0.313996 0.0316003 +214.134 0.311427 0.031822 +214.513 0.309374 0.0277796 +214.891 0.310953 0.02845 +215.27 0.314698 0.0297957 +215.649 0.313525 0.0356094 +216.028 0.310479 0.0333595 +216.406 0.317545 0.0304781 +216.785 0.311719 0.0363158 +217.164 0.314738 0.0403688 +217.542 0.308109 0.0358061 +217.921 0.30758 0.0386051 +218.3 0.311623 0.0332577 +218.678 0.305912 0.0314293 +219.057 0.300317 0.031567 +219.435 0.303921 0.0312197 +219.814 0.297787 0.0312733 +220.192 0.294899 0.0375198 +220.57 0.294884 0.035371 +220.949 0.288752 0.0306086 +221.327 0.298086 0.0353437 +221.705 0.295992 0.0349822 +222.084 0.291475 0.0401682 +222.462 0.292417 0.0366013 +222.84 0.299564 0.0317647 +223.218 0.292982 0.0297147 +223.597 0.303344 0.0295109 +223.975 0.300941 0.0287965 +224.353 0.309024 0.0327674 +224.731 0.310884 0.0308288 +225.109 0.308045 0.0351337 +225.487 0.316169 0.0352102 +225.865 0.321569 0.0347796 +226.243 0.323347 0.0380226 +226.621 0.329447 0.0343911 +226.999 0.327513 0.0380206 +227.377 0.336823 0.0381278 +227.755 0.335278 0.0382196 +228.133 0.343526 0.0341398 +228.511 0.343986 0.0296446 +228.889 0.343431 0.0291393 +229.266 0.341377 0.0308074 +229.644 0.345808 0.0302053 +230.022 0.347545 0.0284951 +230.4 0.348054 0.0286774 +230.777 0.348256 0.032855 +231.155 0.34948 0.0320609 +231.533 0.349117 0.0323248 +231.91 0.354845 0.0296446 +232.288 0.356118 0.0303817 +232.665 0.356677 0.0289302 +233.043 0.366704 0.027364 +233.42 0.366015 0.0263471 +233.798 0.372568 0.0262716 +234.175 0.382063 0.0274226 +234.553 0.387418 0.0270331 +234.93 0.394877 0.0299649 +235.307 0.40475 0.033766 +235.685 0.409095 0.0339657 +236.062 0.415169 0.0335272 +236.439 0.421273 0.0349371 +236.816 0.429126 0.035306 +237.194 0.437766 0.0345561 +237.571 0.438938 0.0332069 +237.948 0.446703 0.0322202 +238.325 0.447531 0.0325013 +238.702 0.454496 0.0308412 +239.079 0.455475 0.03055 +239.456 0.455388 0.0283888 +239.833 0.457917 0.0275723 +240.21 0.465252 0.0270862 +240.587 0.464392 0.026979 +240.964 0.460486 0.0263623 +241.341 0.460277 0.0283078 +241.718 0.459022 0.0295044 +242.095 0.458616 0.0282639 +242.472 0.458388 0.0299158 +242.849 0.457688 0.0284324 +243.225 0.453717 0.028117 +243.602 0.457933 0.027677 +243.979 0.461386 0.0254614 +244.356 0.461422 0.0255465 +244.732 0.461721 0.0257173 +245.109 0.463904 0.0264017 +245.486 0.469576 0.0268126 +245.862 0.473124 0.0286634 +246.239 0.475876 0.0289264 +246.615 0.480743 0.0312862 +246.992 0.48651 0.031287 +247.368 0.496428 0.0333811 +247.745 0.495659 0.0324508 +248.121 0.503605 0.0341907 +248.497 0.508855 0.0327337 +248.874 0.515152 0.0319791 +249.25 0.516334 0.034346 +249.626 0.521693 0.0331971 +250.003 0.52526 0.030371 +250.379 0.525559 0.0308198 +250.755 0.534607 0.0296429 +251.131 0.534446 0.0293879 +251.507 0.535046 0.0280511 +251.884 0.535992 0.0271496 +252.26 0.537248 0.0268609 +252.636 0.533321 0.0266093 +253.012 0.535612 0.0266616 +253.388 0.535937 0.026279 +253.764 0.531982 0.0270695 +254.14 0.528369 0.0270749 +254.516 0.527393 0.0276897 +254.892 0.521728 0.0280492 +255.268 0.523224 0.0275182 +255.643 0.516338 0.0266803 +256.019 0.519041 0.0261889 +256.395 0.512929 0.0258839 +256.771 0.513073 0.0254951 +257.147 0.510577 0.0256262 +257.522 0.512827 0.0257806 +257.898 0.513743 0.0263631 +258.274 0.508347 0.0270949 +258.649 0.514401 0.0279686 +259.025 0.512006 0.0290798 +259.4 0.512755 0.0296772 +259.776 0.514811 0.0307294 +260.151 0.525315 0.0302717 +260.527 0.518848 0.0316213 +260.902 0.526314 0.0324939 +261.278 0.52454 0.0317627 +261.653 0.531422 0.0332104 +262.028 0.533974 0.0334216 +262.404 0.536043 0.0309709 +262.779 0.543352 0.0313819 +263.154 0.54794 0.0307154 +263.53 0.550762 0.0320469 +263.905 0.552528 0.0307398 +264.28 0.554879 0.0295488 +264.655 0.559032 0.0281059 +265.03 0.55911 0.0286681 +265.405 0.566009 0.0274206 +265.781 0.569792 0.0276956 +266.156 0.570956 0.0274132 +266.531 0.567272 0.0277839 +266.906 0.570283 0.0274729 +267.281 0.570774 0.0272539 +267.655 0.569833 0.028007 +268.03 0.568294 0.0278111 +268.405 0.564687 0.0271241 +268.78 0.564949 0.0269975 +269.155 0.562654 0.0261203 +269.53 0.561478 0.027563 +269.904 0.558484 0.0266735 +270.279 0.55748 0.0257763 +270.654 0.550065 0.0256621 +271.029 0.549211 0.0254734 +271.403 0.544193 0.0257006 +271.778 0.545811 0.0258722 +272.152 0.542185 0.0256672 +272.527 0.539289 0.0271266 +272.901 0.540542 0.0286853 +273.276 0.536968 0.0278663 +273.65 0.539334 0.0292708 +274.025 0.530968 0.0297746 +274.399 0.533161 0.0311677 +274.774 0.534784 0.0320431 +275.148 0.53533 0.0325169 +275.522 0.537426 0.0330956 +275.897 0.537044 0.0338911 +276.271 0.53997 0.0302203 +276.645 0.540722 0.033205 +277.019 0.543189 0.0333782 +277.393 0.541555 0.0324658 +277.768 0.547453 0.032261 +278.142 0.551185 0.031921 +278.516 0.555177 0.0316211 +278.89 0.559372 0.0321736 +279.264 0.565597 0.0300815 +279.638 0.566393 0.0308727 +280.012 0.570083 0.0309665 +280.386 0.576132 0.0305337 +280.76 0.577174 0.0287674 +281.133 0.586509 0.0298414 +281.507 0.586814 0.0291817 +281.881 0.588273 0.0285092 +282.255 0.591796 0.0281795 +282.629 0.592962 0.0286035 +283.002 0.592806 0.0291989 +283.376 0.593257 0.0276655 +283.75 0.596407 0.0283997 +284.123 0.596591 0.0286211 +284.497 0.594255 0.0277782 +284.87 0.597469 0.0277295 +285.244 0.59654 0.0290019 +285.618 0.597337 0.0283355 +285.991 0.594493 0.0281928 +286.364 0.596671 0.027795 +286.738 0.592585 0.0267845 +287.111 0.594068 0.0261081 +287.485 0.590085 0.0256243 +287.858 0.591912 0.0254777 +288.231 0.583713 0.0255006 +288.605 0.586952 0.0255562 +288.978 0.582371 0.0257394 +289.351 0.580222 0.0264698 +289.724 0.57348 0.0275021 +290.097 0.576042 0.0271031 +290.47 0.575124 0.028561 +290.843 0.573819 0.0294149 +291.217 0.56926 0.0296457 +291.59 0.567571 0.0297788 +291.963 0.568269 0.0314439 +292.336 0.564272 0.0295361 +292.708 0.566189 0.0316097 +293.081 0.568004 0.0320434 +293.454 0.566287 0.0325621 +293.827 0.565044 0.032524 +294.2 0.56726 0.0315215 +294.573 0.563374 0.0322335 +294.945 0.567884 0.0315973 +295.318 0.568415 0.032445 +295.691 0.573258 0.0327578 +296.064 0.574821 0.0325531 +296.436 0.575293 0.030647 +296.809 0.572967 0.0310556 +297.181 0.576155 0.0305103 +297.554 0.580017 0.0302667 +297.926 0.5846 0.0303502 +298.299 0.58696 0.0289366 +298.671 0.588082 0.0290244 +299.044 0.59683 0.0282693 +299.416 0.597222 0.0290704 +299.788 0.600656 0.0282716 +300.161 0.607342 0.0280118 +300.533 0.610116 0.0283949 +300.905 0.609665 0.0292198 +301.278 0.612309 0.0279523 +301.65 0.615223 0.0274741 +302.022 0.617192 0.0279806 +302.394 0.616474 0.028122 +302.766 0.620764 0.0283016 +303.138 0.622911 0.028303 +303.51 0.625437 0.0286314 +303.882 0.624028 0.0281742 +304.254 0.628581 0.0279468 +304.626 0.628709 0.0280837 +304.998 0.630648 0.0272891 +305.37 0.62646 0.0269211 +305.742 0.627503 0.0275406 +306.114 0.627618 0.0266191 +306.485 0.62421 0.026308 +306.857 0.626586 0.0259915 +307.229 0.626382 0.0256064 +307.601 0.623472 0.0259351 +307.972 0.620705 0.025571 +308.344 0.616586 0.0260381 +308.716 0.61462 0.0258908 +309.087 0.616022 0.0266871 +309.459 0.612635 0.0266319 +309.83 0.614311 0.0268121 +310.202 0.608867 0.0289386 +310.573 0.609149 0.0289214 +310.945 0.607125 0.0301391 +311.316 0.602231 0.0311533 +311.687 0.599452 0.030487 +312.059 0.598129 0.0316568 +312.43 0.595141 0.0312781 +312.801 0.597139 0.0311007 +313.172 0.591056 0.03208 +313.544 0.59224 0.0329595 +313.915 0.584079 0.0345785 +314.286 0.587914 0.0323044 +314.657 0.581095 0.0320146 +315.028 0.582874 0.0335523 +315.399 0.578894 0.0319013 +315.77 0.58159 0.0320141 +316.141 0.578159 0.0314363 +316.512 0.579862 0.0310178 +316.883 0.579401 0.0331201 +317.254 0.577843 0.0330037 +317.625 0.581846 0.0306493 +317.996 0.577482 0.0306178 +318.366 0.578646 0.0306232 +318.737 0.581538 0.029072 +319.108 0.584022 0.0293974 +319.478 0.58249 0.0294665 +319.849 0.586233 0.0304178 +320.22 0.590619 0.0302139 +320.59 0.587575 0.0288091 +320.961 0.590521 0.0287788 +321.331 0.593471 0.0286253 +321.702 0.59444 0.0281329 +322.072 0.601174 0.0278274 +322.443 0.60048 0.0292316 +322.813 0.601782 0.0281045 +323.184 0.605895 0.0289469 +323.554 0.609532 0.0299816 +323.924 0.608847 0.0292753 +324.295 0.616473 0.030602 +324.665 0.61305 0.0285739 +325.035 0.616081 0.0293887 +325.405 0.616758 0.0294862 +325.775 0.628098 0.030477 +326.145 0.619578 0.0308206 +326.515 0.621216 0.029772 +326.886 0.627738 0.0295914 +327.256 0.626568 0.0295658 +327.626 0.629406 0.0296142 +327.995 0.62376 0.0294958 +328.365 0.615375 0.0287549 +328.735 0.611556 0.0288387 +329.105 0.611297 0.0282967 +329.475 0.627015 0.0275947 +329.845 0.649948 0.0262738 +330.215 0.670846 0.025979 +330.584 0.682238 0.0261857 +330.954 0.676892 0.0262245 +331.324 0.657654 0.0257249 +331.693 0.644622 0.0260649 +332.063 0.630791 0.0257951 +332.432 0.618079 0.0257094 +332.802 0.614814 0.0257108 +333.171 0.617131 0.0276727 +333.541 0.617764 0.0289819 +333.91 0.617509 0.0286637 +334.28 0.619797 0.0284351 +334.649 0.61988 0.0301548 +335.018 0.618932 0.0302053 +335.388 0.621832 0.0307945 +335.757 0.617712 0.0295769 +336.126 0.614386 0.0315951 +336.495 0.613389 0.0309798 +336.865 0.610337 0.0318954 +337.234 0.61166 0.0327048 +337.603 0.610667 0.0309621 +337.972 0.60876 0.0327331 +338.341 0.604757 0.0339816 +338.71 0.600771 0.0338327 +339.079 0.602473 0.0327258 +339.448 0.601701 0.032525 +339.817 0.598736 0.0336063 +340.186 0.596368 0.0326231 +340.554 0.593179 0.0335713 +340.923 0.589756 0.0329337 +341.292 0.590175 0.0335235 +341.661 0.590335 0.0321355 +342.029 0.585032 0.0311129 +342.398 0.582937 0.0306218 +342.767 0.583952 0.0313546 +343.135 0.584563 0.0308058 +343.504 0.58375 0.0303267 +343.872 0.583772 0.031949 +344.241 0.585973 0.0319504 +344.609 0.584099 0.0307217 +344.978 0.582581 0.0295322 +345.346 0.583963 0.0303576 +345.715 0.58216 0.029615 +346.083 0.581764 0.0289383 +346.451 0.584829 0.0288906 +346.82 0.583778 0.0289664 +347.188 0.582083 0.028232 +347.556 0.587194 0.0273812 +347.924 0.588094 0.0281717 +348.292 0.589288 0.0279594 +348.66 0.590008 0.0285207 +349.028 0.591043 0.0280383 +349.396 0.594539 0.0289239 +349.764 0.600208 0.0284423 +350.132 0.59376 0.0279001 +350.5 0.598496 0.0295638 +350.868 0.600639 0.0297339 +351.236 0.598262 0.0311967 +351.604 0.602866 0.0301324 +351.972 0.60887 0.0301506 +352.34 0.605464 0.0303712 +352.707 0.610508 0.0303362 +353.075 0.612801 0.0305782 +353.443 0.614758 0.0297859 +353.81 0.613761 0.02915 +354.178 0.620059 0.0291637 +354.545 0.61874 0.0289472 +354.913 0.621893 0.0293033 +355.28 0.618467 0.0286873 +355.648 0.621781 0.0296367 +356.015 0.623872 0.029588 +356.383 0.622358 0.0286674 +356.75 0.624138 0.0277521 +357.117 0.622129 0.0270451 +357.485 0.625561 0.0267695 +357.852 0.631763 0.0265928 +358.219 0.625642 0.0263338 +358.586 0.629777 0.0257707 +358.954 0.629975 0.0256186 +359.321 0.625092 0.025584 +359.688 0.630907 0.025711 +360.055 0.626927 0.0257248 +360.422 0.629591 0.0260143 +360.789 0.629552 0.0260533 +361.156 0.626909 0.0268567 +361.523 0.62527 0.0272637 +361.89 0.626782 0.0287785 +362.256 0.624848 0.0270521 +362.623 0.628796 0.0281498 +362.99 0.622848 0.0294971 +363.357 0.627365 0.0283775 +363.724 0.624341 0.0294912 +364.09 0.617742 0.0305794 +364.457 0.62705 0.0301697 +364.823 0.61967 0.0302328 +365.19 0.614386 0.0304079 +365.557 0.613178 0.0314806 +365.923 0.615492 0.0307249 +366.29 0.618271 0.0312387 +366.656 0.611956 0.0313234 +367.022 0.61317 0.031758 +367.389 0.609597 0.0313103 +367.755 0.608931 0.0325569 +368.121 0.604419 0.0329693 +368.488 0.599927 0.0333944 +368.854 0.602051 0.0334422 +369.22 0.602837 0.0329387 +369.586 0.595205 0.0329482 +369.952 0.595254 0.0327899 +370.318 0.590219 0.0322267 +370.684 0.590757 0.0323716 +371.051 0.589849 0.0328238 +371.416 0.590177 0.0326765 +371.782 0.584478 0.0322089 +372.148 0.585202 0.0321296 +372.514 0.584309 0.0324137 +372.88 0.580982 0.0319921 +373.246 0.581086 0.0319318 +373.612 0.579962 0.0300687 +373.977 0.58118 0.0299098 +374.343 0.575422 0.0297319 +374.709 0.580397 0.0288123 +375.074 0.573142 0.0297127 +375.44 0.572463 0.0303409 +375.806 0.574874 0.0297494 +376.171 0.568458 0.0293494 +376.537 0.561565 0.0284781 +376.902 0.565599 0.029271 +377.268 0.569739 0.0284894 +377.633 0.568274 0.0285874 +377.998 0.565066 0.0282925 +378.364 0.571603 0.0281253 +378.729 0.568462 0.0279737 +379.094 0.575998 0.027391 +379.459 0.578966 0.0278692 +379.825 0.568078 0.0272889 +380.19 0.569828 0.028339 +380.555 0.569781 0.0275321 +380.92 0.576475 0.0280239 +381.285 0.574087 0.0293121 +381.65 0.568912 0.0291651 +382.015 0.571627 0.0291197 +382.38 0.573685 0.0286327 +382.745 0.577935 0.0292801 +383.11 0.579401 0.0286528 +383.474 0.578578 0.0302708 +383.839 0.580786 0.0291325 +384.204 0.585235 0.0302693 +384.569 0.584739 0.0299607 +384.933 0.591544 0.0311716 +385.298 0.588022 0.0311823 +385.663 0.583911 0.0298461 +386.027 0.585716 0.0307457 +386.392 0.58565 0.0307117 +386.756 0.585616 0.0297866 +387.121 0.59152 0.0286595 +387.485 0.598445 0.0306746 +387.849 0.596746 0.0288751 +388.214 0.603594 0.0286687 +388.578 0.602131 0.0292865 +388.942 0.603489 0.0277324 +389.307 0.598692 0.0286667 +389.671 0.607781 0.0282973 +390.035 0.60522 0.0281641 +390.399 0.605572 0.0274739 +390.763 0.607045 0.0265432 +391.127 0.607395 0.0263146 +391.492 0.608139 0.0258237 +391.856 0.608557 0.0260667 +392.219 0.610198 0.0261217 +392.583 0.617168 0.0256216 +392.947 0.626732 0.0256953 +393.311 0.617092 0.0258901 +393.675 0.617535 0.0260006 +394.039 0.618466 0.0260713 +394.403 0.623267 0.0261008 +394.766 0.619251 0.0269313 +395.13 0.618664 0.0261293 +395.494 0.615242 0.0265709 +395.857 0.616992 0.0276001 +396.221 0.616277 0.0280275 +396.584 0.614027 0.0275111 +396.948 0.612448 0.0289128 +397.311 0.612894 0.0277412 +397.675 0.615562 0.0297034 +398.038 0.618332 0.0291495 +398.401 0.624496 0.0297859 +398.765 0.620363 0.0293702 +399.128 0.625873 0.0297153 +399.491 0.61585 0.0310453 +399.854 0.610532 0.0308681 +400.218 0.617293 0.0322966 +400.581 0.615604 0.0308607 +400.944 0.616536 0.0314934 +401.307 0.620034 0.0307587 +401.67 0.613599 0.0305735 +402.033 0.612888 0.0338462 +402.396 0.613354 0.0305073 +402.759 0.607713 0.0328965 +403.122 0.606443 0.032103 +403.485 0.611229 0.0333828 +403.847 0.601353 0.0328073 +404.21 0.603427 0.0323272 +404.573 0.607499 0.0338893 +404.936 0.604207 0.0319037 +405.298 0.612727 0.0346819 +405.661 0.605865 0.0335537 +406.023 0.606417 0.0327084 +406.386 0.601479 0.0358988 +406.749 0.601842 0.0331432 +407.111 0.59455 0.0325008 +407.473 0.585777 0.0336173 +407.836 0.590621 0.0327326 +408.198 0.590611 0.0333083 +408.561 0.591098 0.0318874 +408.923 0.587236 0.0337458 +409.285 0.584024 0.0328435 +409.647 0.588361 0.0319653 +410.009 0.588694 0.0307217 +410.372 0.582749 0.0308982 +410.734 0.583888 0.0304415 +411.096 0.57805 0.0306533 +411.458 0.574059 0.0292396 +411.82 0.573261 0.0310782 +412.182 0.57118 0.0297501 +412.544 0.57052 0.0299512 +412.906 0.572526 0.0295077 +413.267 0.572396 0.0300861 +413.629 0.572412 0.0286031 +413.991 0.571384 0.0288768 +414.353 0.569729 0.0279909 +414.714 0.569042 0.0288146 +415.076 0.561213 0.0292568 +415.438 0.565865 0.0279677 +415.799 0.561748 0.028221 +416.161 0.562588 0.0283598 +416.522 0.560522 0.0280512 +416.884 0.556682 0.0277771 +417.245 0.564384 0.0269462 +417.607 0.556233 0.0277111 +417.968 0.549893 0.0280484 +418.329 0.557883 0.0275662 +418.691 0.552522 0.0282038 +419.052 0.553893 0.0283656 +419.413 0.544975 0.0283944 +419.774 0.55274 0.0287474 +420.135 0.544631 0.0294337 +420.497 0.544025 0.0283494 +420.858 0.544368 0.0283741 +421.219 0.556108 0.0294343 +421.58 0.558275 0.0284173 +421.941 0.563275 0.0292404 +422.302 0.556568 0.0305945 +422.662 0.558035 0.0307282 +423.023 0.558024 0.0297561 +423.384 0.556993 0.0303786 +423.745 0.554062 0.0306244 +424.106 0.558742 0.0307116 +424.466 0.559095 0.0312809 +424.827 0.568302 0.0301537 +425.187 0.568082 0.0307669 +425.548 0.561527 0.031037 +425.909 0.562144 0.0298568 +426.269 0.56049 0.0299316 +426.63 0.558856 0.0318824 +426.99 0.562572 0.0303461 +427.35 0.567362 0.0289206 +427.711 0.561463 0.0286239 +428.071 0.561738 0.0291759 +428.431 0.554199 0.0296167 +428.792 0.555323 0.0287825 +429.152 0.561977 0.0289679 +429.512 0.571685 0.029322 +429.872 0.573348 0.0280719 +430.232 0.582303 0.0274407 +430.592 0.592775 0.0266035 +430.952 0.590523 0.0266319 +431.312 0.59495 0.027305 +431.672 0.578086 0.0273391 +432.032 0.58229 0.0266664 +432.392 0.575612 0.026602 +432.752 0.56076 0.026213 +433.111 0.573482 0.0265922 +433.471 0.598149 0.025874 +433.831 0.608249 0.0260217 +434.19 0.591065 0.0258724 +434.55 0.589087 0.026188 +434.91 0.589976 0.0260406 +435.269 0.586867 0.025913 +435.629 0.58812 0.025897 +435.988 0.590619 0.0259106 +436.348 0.600831 0.0258956 +436.707 0.601808 0.0260627 +437.066 0.592393 0.0266962 +437.426 0.587905 0.0262811 +437.785 0.590485 0.0264292 +438.144 0.594598 0.0275174 +438.503 0.602167 0.0268779 +438.862 0.602189 0.027512 +439.222 0.604232 0.0285958 +439.581 0.600449 0.0278131 +439.94 0.590429 0.0283447 +440.299 0.585891 0.0279933 +440.658 0.587507 0.0279212 +441.017 0.603078 0.0285702 +441.376 0.616915 0.0288598 +441.734 0.609922 0.0290985 +442.093 0.60062 0.0308725 +442.452 0.592893 0.0302365 +442.811 0.600934 0.029905 +443.169 0.603393 0.0322014 +443.528 0.613656 0.0306168 +443.887 0.603891 0.0322245 +444.245 0.60646 0.0315327 +444.604 0.604727 0.0318874 +444.962 0.599031 0.0324345 +445.321 0.599896 0.0308681 +445.679 0.594657 0.032113 +446.038 0.602737 0.0313387 +446.396 0.609336 0.0309624 +446.754 0.606781 0.0338131 +447.112 0.608544 0.0341903 +447.471 0.599768 0.0326096 +447.829 0.6012 0.0348544 +448.187 0.600359 0.0317914 +448.545 0.59918 0.0336286 +448.903 0.595149 0.0326638 +449.261 0.601911 0.0338099 +449.619 0.598425 0.0326641 +449.977 0.603765 0.0354624 +450.335 0.606897 0.0366003 +450.693 0.602507 0.0345538 +451.051 0.592405 0.0354352 +451.409 0.595151 0.034062 +451.766 0.600877 0.03291 +452.124 0.603029 0.0338455 +452.482 0.598956 0.0344594 +452.839 0.599459 0.0341529 +453.197 0.592825 0.0336488 +453.555 0.595983 0.0345621 +453.912 0.591542 0.0346212 +454.27 0.579561 0.03353 +454.627 0.586237 0.0345358 +454.984 0.59018 0.0325867 +455.342 0.593235 0.0341864 +455.699 0.593847 0.032446 +456.056 0.580235 0.0326976 +456.414 0.58312 0.030559 +456.771 0.569605 0.0329904 +457.128 0.580349 0.0294419 +457.485 0.575538 0.0304396 +457.842 0.583304 0.0306437 +458.199 0.580716 0.0313972 +458.556 0.577178 0.0311521 +458.913 0.576314 0.0307559 +459.27 0.576849 0.0304996 +459.627 0.577391 0.029021 +459.984 0.570078 0.0293733 +460.341 0.572136 0.0299956 +460.697 0.571466 0.0287204 +461.054 0.582195 0.0284166 +461.411 0.57319 0.0290441 +461.767 0.571563 0.027762 +462.124 0.556345 0.0287709 +462.481 0.5589 0.0285597 +462.837 0.55391 0.0292585 +463.194 0.554463 0.0280362 +463.55 0.554293 0.0275045 +463.907 0.547446 0.0281903 +464.263 0.548333 0.0267337 +464.619 0.550594 0.0271457 +464.976 0.548883 0.0288423 +465.332 0.558864 0.0278493 +465.688 0.57298 0.0272285 +466.044 0.569605 0.0287676 +466.4 0.559746 0.0285885 +466.756 0.551938 0.0288606 +467.112 0.55093 0.027745 +467.468 0.54938 0.0286056 +467.824 0.556859 0.0286624 +468.18 0.552141 0.028546 +468.536 0.543579 0.0280772 +468.892 0.549875 0.0280809 +469.248 0.547783 0.0287935 +469.604 0.544863 0.0292464 +469.959 0.547348 0.029703 +470.315 0.545215 0.0300781 +470.671 0.546533 0.0294304 +471.026 0.53641 0.0300466 +471.382 0.532465 0.0302432 +471.737 0.535027 0.0299159 +472.093 0.533828 0.0300216 +472.448 0.528625 0.0298291 +472.804 0.530491 0.0306727 +473.159 0.534323 0.0309984 +473.514 0.541829 0.0310279 +473.87 0.531457 0.0308895 +474.225 0.525917 0.0339813 +474.58 0.52828 0.0309263 +474.935 0.543661 0.0294789 +475.29 0.549873 0.0306252 +475.645 0.537851 0.0298794 +476 0.533686 0.0313185 +476.355 0.54266 0.0324496 +476.71 0.5354 0.0294599 +477.065 0.521153 0.030738 +477.42 0.527969 0.031286 +477.775 0.535449 0.0309279 +478.13 0.543473 0.0299149 +478.485 0.539155 0.029577 +478.839 0.528792 0.0325537 +479.194 0.527202 0.0287109 +479.549 0.532539 0.0303795 +479.903 0.54718 0.0293425 +480.258 0.54314 0.0287952 +480.612 0.538432 0.0292989 +480.967 0.532301 0.0280063 +481.321 0.542812 0.0274369 +481.675 0.553183 0.0270433 +482.03 0.552884 0.0271109 +482.384 0.541853 0.0269582 +482.738 0.538107 0.0273128 +483.093 0.539858 0.0271405 +483.447 0.53746 0.0269201 +483.801 0.531525 0.026897 +484.155 0.54169 0.0259622 +484.509 0.457682 0.027123 +484.863 0.434371 0.0286092 +485.217 0.511776 0.026942 +485.571 0.600074 0.0254385 +485.925 0.709644 0.0264408 +486.279 0.636706 0.0256999 +486.632 0.56873 0.0257318 +486.986 0.568745 0.0255715 +487.34 0.579694 0.0260481 +487.694 0.569205 0.0259397 +488.047 0.555782 0.0259178 +488.401 0.563499 0.0259622 +488.754 0.565026 0.0263402 +489.108 0.569376 0.0262615 +489.461 0.562981 0.0257855 +489.815 0.571585 0.0262361 +490.168 0.560832 0.0257546 +490.522 0.565318 0.0259031 +490.875 0.567107 0.0264482 +491.228 0.576435 0.0271279 +491.581 0.570987 0.0264324 +491.935 0.575453 0.0269476 +492.288 0.57152 0.026694 +492.641 0.573311 0.0261918 +492.994 0.570534 0.0263629 +493.347 0.571685 0.0264431 +493.7 0.574425 0.0269525 +494.053 0.575266 0.0274378 +494.406 0.586956 0.0279797 +494.759 0.578531 0.029164 +495.111 0.560688 0.0282904 +495.464 0.558005 0.0278684 +495.817 0.57059 0.0276233 +496.169 0.596183 0.0292932 +496.522 0.604252 0.028846 +496.875 0.59269 0.0291011 +497.227 0.581562 0.0291912 +497.58 0.573125 0.0304647 +497.932 0.576792 0.0294988 +498.285 0.575108 0.0298959 +498.637 0.570421 0.0298186 +498.989 0.576243 0.0311536 +499.342 0.588289 0.0300187 +499.694 0.597422 0.0310836 +500.046 0.58384 0.0301021 +500.398 0.57326 0.0313043 +500.751 0.586768 0.030633 +501.103 0.583617 0.030465 +501.455 0.588139 0.0315444 +501.807 0.587685 0.0307228 +502.159 0.586555 0.0308741 +502.511 0.587915 0.0316328 +502.862 0.585961 0.0330016 +503.214 0.590471 0.0340922 +503.566 0.595161 0.0325143 +503.918 0.585345 0.0347426 +504.27 0.583884 0.0335272 +504.621 0.587913 0.0325484 +504.973 0.592267 0.0337636 +505.324 0.587229 0.0341225 +505.676 0.57893 0.032549 +506.027 0.587287 0.0337624 +506.379 0.600027 0.034581 +506.73 0.586466 0.03298 +507.082 0.579162 0.0341952 +507.433 0.574066 0.0331743 +507.784 0.588473 0.0364293 +508.136 0.588946 0.0342998 +508.487 0.593432 0.0331521 +508.838 0.584264 0.0354166 +509.189 0.582723 0.0356851 +509.54 0.57723 0.0319887 +509.891 0.587124 0.0352833 +510.242 0.584271 0.0330177 +510.593 0.584658 0.0328637 +510.944 0.581464 0.0344749 +511.295 0.591176 0.0333661 +511.646 0.59155 0.0324034 +511.996 0.573702 0.0335215 +512.347 0.57538 0.0331819 +512.698 0.587017 0.0327461 +513.049 0.590724 0.0308434 +513.399 0.585615 0.0345095 +513.75 0.578997 0.0331207 +514.1 0.591101 0.0316035 +514.451 0.592594 0.0307238 +514.801 0.574881 0.0317606 +515.152 0.574815 0.0301753 +515.502 0.572352 0.0319868 +515.852 0.575777 0.0312308 +516.202 0.57891 0.0297258 +516.553 0.591031 0.0310853 +516.903 0.582046 0.0303495 +517.253 0.58069 0.0313237 +517.603 0.56886 0.0303966 +517.953 0.570996 0.0304227 +518.303 0.565432 0.0287098 +518.653 0.569803 0.0287139 +519.003 0.58186 0.0293606 +519.353 0.586945 0.0305296 +519.703 0.577304 0.0298749 +520.052 0.57067 0.028465 +520.402 0.565291 0.0299911 +520.752 0.570587 0.0296169 +521.102 0.580395 0.0286987 +521.451 0.578073 0.0304984 +521.801 0.568626 0.0285153 +522.15 0.559511 0.0284314 +522.5 0.572333 0.0266584 +522.849 0.563608 0.0286534 +523.199 0.563974 0.0272458 +523.548 0.55834 0.0282131 +523.897 0.566817 0.0277995 +524.246 0.571761 0.0273121 +524.596 0.565857 0.0287429 +524.945 0.55533 0.0266931 +525.294 0.558448 0.02826 +525.643 0.556705 0.0278293 +525.992 0.558806 0.0270698 +526.341 0.555636 0.0284932 +526.69 0.557286 0.0280306 +527.039 0.555645 0.0277115 +527.388 0.548432 0.0272612 +527.737 0.53781 0.0277962 +528.086 0.534314 0.0294656 +528.434 0.550318 0.0284053 +528.783 0.567602 0.0281956 +529.132 0.5576 0.0283723 +529.48 0.548741 0.02851 +529.829 0.549848 0.0279172 +530.177 0.553973 0.0279917 +530.526 0.544499 0.0276912 +530.874 0.54175 0.0308272 +531.223 0.538539 0.0293261 +531.571 0.545343 0.0292626 +531.919 0.536838 0.0291813 +532.268 0.520942 0.0290984 +532.616 0.523235 0.0309625 +532.964 0.537869 0.0302568 +533.312 0.538573 0.0287143 +533.66 0.537389 0.0310693 +534.008 0.5341 0.0305404 +534.356 0.534345 0.0303121 +534.704 0.540351 0.0294202 +535.052 0.540637 0.0297397 +535.4 0.538095 0.0308079 +535.748 0.534682 0.0288309 +536.095 0.533983 0.0291252 +536.443 0.523539 0.0304092 +536.791 0.513707 0.0303444 +537.139 0.516464 0.0309812 +537.486 0.518218 0.0302636 +537.834 0.522192 0.031066 +538.181 0.52488 0.0298733 +538.529 0.527447 0.031819 +538.876 0.52093 0.0321387 +539.223 0.528463 0.0287564 +539.571 0.530978 0.0303898 +539.918 0.519689 0.0312235 +540.265 0.5129 0.0306958 +540.613 0.516611 0.0320813 +540.96 0.525154 0.030204 +541.307 0.514041 0.0305069 +541.654 0.515215 0.030393 +542.001 0.517106 0.0297486 +542.348 0.517031 0.0295823 +542.695 0.519881 0.0304871 +543.042 0.516203 0.0289142 +543.388 0.511914 0.0305278 +543.735 0.518624 0.0286599 +544.082 0.511865 0.0304701 +544.429 0.517243 0.029338 +544.775 0.521736 0.0281933 +545.122 0.528067 0.0281698 +545.469 0.519783 0.0296141 +545.815 0.510439 0.0287709 +546.162 0.509308 0.0277355 +546.508 0.507784 0.0284281 +546.854 0.503069 0.0268259 +547.201 0.502074 0.0285578 +547.547 0.525265 0.0285329 +547.893 0.527484 0.0274601 +548.24 0.542698 0.0276487 +548.586 0.529197 0.0271092 +548.932 0.516427 0.0273345 +549.278 0.507789 0.0274756 +549.624 0.516459 0.0277208 +549.97 0.523619 0.0277542 +550.316 0.524349 0.0281139 +550.662 0.523027 0.0264288 +551.008 0.495151 0.0269624 +551.353 0.505734 0.0268523 +551.699 0.519199 0.027043 +552.045 0.541672 0.0258912 +552.391 0.529612 0.0269965 +552.736 0.51963 0.0261153 +553.082 0.512515 0.0265241 +553.427 0.532739 0.0262296 +553.773 0.531855 0.0267942 +554.118 0.517272 0.0264223 +554.464 0.508571 0.0268976 +554.809 0.506885 0.0263006 +555.154 0.524313 0.0260774 +555.5 0.532519 0.0264412 +555.845 0.541142 0.0263153 +556.19 0.533066 0.0271187 +556.535 0.530279 0.0264977 +556.88 0.525004 0.0261448 +557.225 0.532907 0.0261406 +557.57 0.54101 0.0267004 +557.915 0.542466 0.0270221 +558.26 0.5388 0.0261761 +558.605 0.532268 0.0263283 +558.95 0.525364 0.0265716 +559.295 0.517689 0.0261204 +559.639 0.526582 0.0268433 +559.984 0.525948 0.0262388 +560.329 0.538263 0.0261058 +560.673 0.540503 0.0271135 +561.018 0.547785 0.0268558 +561.362 0.541118 0.0263366 +561.707 0.542807 0.0265545 +562.051 0.535206 0.0271576 +562.395 0.52916 0.0268643 +562.74 0.538148 0.0264809 +563.084 0.540013 0.026696 +563.428 0.544777 0.0272888 +563.772 0.551159 0.0265905 +564.117 0.541434 0.0269706 +564.461 0.545237 0.0273358 +564.805 0.540771 0.0272148 +565.149 0.542532 0.0281509 +565.493 0.547976 0.0286647 +565.836 0.537211 0.0276164 +566.18 0.541974 0.028391 +566.524 0.526083 0.0265081 +566.868 0.547591 0.0277788 +567.212 0.552399 0.0276411 +567.555 0.544718 0.0284684 +567.899 0.534647 0.0282838 +568.242 0.539 0.0284443 +568.586 0.533035 0.0281112 +568.929 0.531913 0.0277405 +569.273 0.543805 0.028527 +569.616 0.55871 0.0292119 +569.96 0.545193 0.028707 +570.303 0.543001 0.0288128 +570.646 0.547149 0.029559 +570.989 0.540681 0.0293023 +571.333 0.56054 0.0294357 +571.676 0.568568 0.0299419 +572.019 0.564463 0.0309047 +572.362 0.521707 0.0292217 +572.705 0.527244 0.0297027 +573.048 0.551126 0.0303122 +573.391 0.550749 0.0301171 +573.733 0.555404 0.0313933 +574.076 0.557682 0.031012 +574.419 0.57193 0.0316461 +574.762 0.56388 0.0322438 +575.104 0.570039 0.0321622 +575.447 0.531032 0.0314852 +575.789 0.530018 0.0315881 +576.132 0.553112 0.0313287 +576.474 0.551706 0.0326602 +576.817 0.543192 0.0309099 +577.159 0.558996 0.030796 +577.501 0.593744 0.0316625 +577.844 0.593891 0.0327423 +578.186 0.576874 0.0337802 +578.528 0.543078 0.0340429 +578.87 0.546166 0.033107 +579.212 0.53801 0.0324341 +579.554 0.561401 0.0331082 +579.896 0.561847 0.0341942 +580.238 0.533709 0.0333043 +580.58 0.549075 0.0320922 +580.922 0.600255 0.0322099 +581.264 0.648919 0.0326095 +581.606 0.605679 0.0331188 +581.947 0.58514 0.0341941 +582.289 0.575336 0.0340567 +582.631 0.560765 0.0336747 +582.972 0.545368 0.0335865 +583.314 0.542034 0.0335234 +583.655 0.538566 0.0329895 +583.997 0.537869 0.0327311 +584.338 0.555248 0.033224 +584.679 0.569105 0.0339096 +585.021 0.592487 0.0339083 +585.362 0.578177 0.0344183 +585.703 0.570558 0.0333564 +586.044 0.532656 0.0333607 +586.385 0.540559 0.0328215 +586.726 0.557876 0.0335218 +587.067 0.589374 0.0333249 +587.408 0.574127 0.0344377 +587.749 0.570218 0.0333658 +588.09 0.564662 0.0345071 +588.431 0.573744 0.0336583 +588.772 0.556308 0.0333823 +589.113 0.550719 0.0332589 +589.453 0.562506 0.0320864 +589.794 0.577414 0.031746 +590.134 0.573368 0.033314 +590.475 0.562502 0.0326102 +590.815 0.55326 0.0323493 +591.156 0.569357 0.0323281 +591.496 0.538795 0.0317864 +591.837 0.532307 0.0303039 +592.177 0.548748 0.0317476 +592.517 0.577646 0.0315384 +592.857 0.589947 0.0331455 +593.198 0.582089 0.0322622 +593.538 0.57165 0.032219 +593.878 0.548205 0.0313963 +594.218 0.536105 0.0305354 +594.558 0.558081 0.0304857 +594.898 0.578679 0.030413 +595.237 0.583187 0.0314057 +595.577 0.568658 0.0299942 +595.917 0.5516 0.0300313 +596.257 0.543098 0.0303745 +596.596 0.546718 0.0296606 +596.936 0.574906 0.0301345 +597.276 0.569621 0.0301426 +597.615 0.545813 0.0279321 +597.955 0.524779 0.029276 +598.294 0.560423 0.0284498 +598.633 0.580512 0.0290981 +598.973 0.568481 0.029541 +599.312 0.546734 0.0277377 +599.651 0.565225 0.0278999 +599.991 0.56136 0.0289123 +600.33 0.56003 0.0283957 +600.669 0.558132 0.0278976 +601.008 0.566852 0.028144 +601.347 0.550628 0.0275752 +601.686 0.545268 0.0277104 +602.025 0.554203 0.0272107 +602.364 0.544612 0.0274413 +602.702 0.55406 0.0275228 +603.041 0.567759 0.0264673 +603.38 0.56598 0.0263241 +603.719 0.562809 0.0269286 +604.057 0.543267 0.0275377 +604.396 0.541082 0.027315 +604.734 0.526525 0.0275347 +605.073 0.551541 0.0293121 +605.411 0.560887 0.0267587 +605.75 0.574163 0.0265468 +606.088 0.550312 0.0267331 +606.426 0.53981 0.0281734 +606.764 0.545505 0.0272206 +607.103 0.565545 0.0270511 +607.441 0.548848 0.0271556 +607.779 0.542819 0.0272731 +608.117 0.541765 0.0276578 +608.455 0.55918 0.0261086 +608.793 0.531337 0.0271033 +609.131 0.521114 0.0300728 +609.469 0.544874 0.0291123 +609.806 0.567052 0.0271996 +610.144 0.557734 0.0264525 +610.482 0.53203 0.0279777 +610.819 0.544443 0.0287266 +611.157 0.543815 0.029117 +611.495 0.543726 0.0266296 +611.832 0.532881 0.0291094 +612.17 0.530057 0.0300007 +612.507 0.540162 0.0299053 +612.844 0.536453 0.0280727 +613.182 0.534554 0.0292558 +613.519 0.533043 0.0294471 +613.856 0.538681 0.0295048 +614.193 0.528286 0.0292612 +614.53 0.529125 0.0297185 +614.867 0.543416 0.0311225 +615.205 0.549993 0.0282433 +615.541 0.54428 0.0290685 +615.878 0.523364 0.0293386 +616.215 0.523079 0.0308009 +616.552 0.528469 0.031073 +616.889 0.524639 0.0325002 +617.226 0.535641 0.0320178 +617.562 0.546515 0.0304188 +617.899 0.54367 0.0296801 +618.235 0.532471 0.0305319 +618.572 0.523496 0.029751 +618.908 0.513997 0.030898 +619.245 0.519961 0.0327324 +619.581 0.530529 0.0322271 +619.918 0.560942 0.0293532 +620.254 0.544661 0.0281792 +620.59 0.506055 0.0306639 +620.926 0.500845 0.0368327 +621.262 0.515821 0.0332694 +621.599 0.540907 0.0316171 +621.935 0.518884 0.0313429 +622.271 0.517471 0.0324037 +622.606 0.520289 0.0321426 +622.942 0.533416 0.0303561 +623.278 0.516876 0.0288668 +623.614 0.502507 0.0345292 +623.95 0.52067 0.0345603 +624.285 0.528383 0.0312082 +624.621 0.525031 0.0295098 +624.957 0.512398 0.0309296 +625.292 0.524617 0.0291898 +625.628 0.513617 0.0310944 +625.963 0.504308 0.0321104 +626.299 0.51126 0.0296732 +626.634 0.515043 0.033852 +626.969 0.524219 0.0303098 +627.304 0.533806 0.0283976 +627.64 0.518199 0.0293404 +627.975 0.49429 0.0320672 +628.31 0.483867 0.031999 +628.645 0.521704 0.0326495 +628.98 0.521587 0.0293103 +629.315 0.518585 0.0295761 +629.65 0.513744 0.0303146 +629.985 0.515391 0.0278644 +630.319 0.518058 0.0291086 +630.654 0.497669 0.0284634 +630.989 0.5075 0.0312091 +631.323 0.519404 0.0285535 +631.658 0.529728 0.03072 +631.993 0.517636 0.0294836 +632.327 0.512526 0.030801 +632.662 0.513321 0.0281132 +632.996 0.505262 0.0293166 +633.33 0.506638 0.0279167 +633.665 0.524265 0.0280619 +633.999 0.517178 0.0283762 +634.333 0.508101 0.0287634 +634.667 0.507856 0.0310666 +635.001 0.510776 0.028769 +635.335 0.50206 0.0304389 +635.669 0.505895 0.0274468 +636.003 0.521541 0.0289799 +636.337 0.527577 0.0288701 +636.671 0.515191 0.0274241 +637.005 0.511817 0.027312 +637.339 0.492841 0.0288439 +637.672 0.497825 0.028606 +638.006 0.510367 0.0293746 +638.339 0.516433 0.0273537 +638.673 0.515398 0.02837 +639.006 0.509903 0.0279578 +639.34 0.504744 0.0283899 +639.673 0.511888 0.0293529 +640.007 0.511546 0.0303158 +640.34 0.506302 0.0280911 +640.673 0.513944 0.0275635 +641.006 0.518369 0.0283797 +641.34 0.518188 0.0283128 +641.673 0.512905 0.0283386 +642.006 0.511675 0.0277363 +642.339 0.515355 0.0289461 +642.672 0.518538 0.0282254 +643.004 0.511298 0.0278876 +643.337 0.530595 0.028205 +643.67 0.511156 0.0278648 +644.003 0.503234 0.027681 +644.336 0.512932 0.0279785 +644.668 0.510381 0.0282127 +645.001 0.535192 0.0287471 +645.333 0.535286 0.0286531 +645.666 0.517353 0.0284838 +645.998 0.520264 0.0289792 +646.331 0.533284 0.0283668 +646.663 0.524578 0.028363 +646.995 0.505662 0.0278043 +647.327 0.507585 0.0276447 +647.66 0.517558 0.0272868 +647.992 0.519197 0.0278636 +648.324 0.512938 0.0280883 +648.656 0.520589 0.0272589 +648.988 0.530194 0.0270199 +649.32 0.540285 0.0275405 +649.652 0.519561 0.0279827 +649.983 0.509915 0.0274962 +650.315 0.52537 0.0279959 +650.647 0.503466 0.0281934 +650.979 0.509809 0.0272202 +651.31 0.505683 0.0275089 +651.642 0.515198 0.0272145 +651.973 0.521666 0.0287276 +652.305 0.517104 0.0275257 +652.636 0.524446 0.0276215 +652.967 0.518129 0.0271044 +653.299 0.516518 0.0286288 +653.63 0.492861 0.0275818 +653.961 0.480062 0.0280266 +654.292 0.402487 0.0282152 +654.623 0.303574 0.0303065 +654.955 0.89983 0.025234 +655.286 0.899775 0.025234 +655.616 0.899719 0.025234 +655.947 0.899662 0.025234 +656.278 0.815899 0.0252123 +656.609 0.345738 0.0439928 +656.94 0.316006 0.0285391 +657.27 0.291108 0.0270551 +657.601 0.300514 0.0268686 +657.932 0.331907 0.0272779 +658.262 0.367305 0.0281031 +658.593 0.417579 0.0275545 +658.923 0.455627 0.0287521 +659.254 0.500217 0.0279855 +659.584 0.52133 0.029625 +659.914 0.534629 0.0285519 +660.244 0.550614 0.0320298 +660.575 0.528537 0.0290046 +660.905 0.544704 0.0307943 +661.235 0.549742 0.0312574 +661.565 0.540902 0.0289022 +661.895 0.55147 0.0284834 +662.225 0.549687 0.0284153 +662.555 0.56179 0.0285186 +662.884 0.536202 0.0294426 +663.214 0.546992 0.0286288 +663.544 0.560903 0.0297768 +663.873 0.550523 0.0298355 +664.203 0.55315 0.0289195 +664.533 0.535818 0.0286876 +664.862 0.560597 0.0290757 +665.192 0.555508 0.0314116 +665.521 0.563987 0.0314261 +665.85 0.539604 0.0305011 +666.18 0.567275 0.0305762 +666.509 0.557871 0.0317 +666.838 0.546625 0.0328848 +667.167 0.562565 0.0312688 +667.496 0.574306 0.0294938 +667.825 0.5722 0.0299231 +668.154 0.572356 0.0296051 +668.483 0.569725 0.0335372 +668.812 0.566784 0.0316438 +669.141 0.57924 0.0347315 +669.47 0.572203 0.0298582 +669.798 0.578758 0.0318092 +670.127 0.555621 0.0330679 +670.456 0.531687 0.030801 +670.784 0.528287 0.0295983 +671.113 0.54216 0.0297424 +671.441 0.563571 0.0295265 +671.77 0.564503 0.0306323 +672.098 0.568496 0.0318393 +672.426 0.584276 0.0365649 +672.755 0.551655 0.0340156 +673.083 0.54432 0.0312202 +673.411 0.545355 0.0316374 +673.739 0.567134 0.0344549 +674.067 0.56209 0.0315475 +674.395 0.562378 0.031256 +674.723 0.557265 0.0296615 +675.051 0.571471 0.0313059 +675.379 0.576414 0.0337818 +675.706 0.590571 0.0326372 +676.034 0.574414 0.0340387 +676.362 0.566568 0.030819 +676.689 0.555748 0.0317172 +677.017 0.566557 0.0304865 +677.344 0.568559 0.0316583 +677.672 0.587676 0.0348072 +677.999 0.558219 0.03447 +678.327 0.556668 0.0303731 +678.654 0.552634 0.0300567 +678.981 0.56036 0.0321983 +679.308 0.577587 0.0313358 +679.636 0.589339 0.0322288 +679.963 0.576144 0.0355391 +680.29 0.54987 0.0348108 +680.617 0.553051 0.0340955 +680.944 0.560806 0.0299524 +681.271 0.593212 0.0360766 +681.597 0.582534 0.0353895 +681.924 0.581389 0.0321303 +682.251 0.55545 0.03545 +682.578 0.554441 0.0324915 +682.904 0.566373 0.0298272 +683.231 0.567215 0.0319277 +683.557 0.5716 0.0348866 +683.884 0.564194 0.0352068 +684.21 0.565826 0.0323891 +684.536 0.558666 0.0316958 +684.863 0.56682 0.0347299 +685.189 0.579019 0.0348519 +685.515 0.575909 0.0334868 +685.841 0.560264 0.0391916 +686.167 0.54042 0.0336864 +686.493 0.572719 0.0321601 +686.819 0.573858 0.0324558 +687.145 0.580903 0.0347647 +687.471 0.579298 0.034018 +687.797 0.574049 0.0343607 +688.123 0.558694 0.0325448 +688.448 0.558241 0.0355142 +688.774 0.577613 0.0337761 +689.1 0.584179 0.0372798 +689.425 0.585839 0.0366122 +689.751 0.574824 0.0354602 +690.076 0.554624 0.0363331 +690.402 0.558762 0.0349427 +690.727 0.576863 0.0311034 +691.052 0.591901 0.0346946 +691.377 0.576204 0.036533 +691.703 0.587318 0.0362756 +692.028 0.566137 0.0343163 +692.353 0.568509 0.0310289 +692.678 0.583909 0.0328081 +693.003 0.594706 0.0326438 +693.328 0.592028 0.0431558 +693.652 0.579223 0.0350719 +693.977 0.552036 0.0410654 +694.302 0.557249 0.0333685 +694.627 0.599095 0.0338924 +694.951 0.602563 0.0381355 +695.276 0.614159 0.0316149 +695.6 0.578695 0.0373285 +695.925 0.53211 0.0337658 +696.249 0.553506 0.0346787 +696.574 0.590679 0.0325109 +696.898 0.613053 0.0346943 +697.222 0.588025 0.0330013 +697.546 0.581246 0.0351651 +697.871 0.567316 0.0316144 +698.195 0.571133 0.0304063 +698.519 0.572404 0.0309585 +698.843 0.595469 0.0321848 +699.167 0.587674 0.0345415 +699.491 0.5862 0.0341568 +699.814 0.578889 0.0327996 +700.138 0.582339 0.0299979 +700.462 0.592669 0.0319626 +700.786 0.586529 0.0305666 +701.109 0.580639 0.0304877 +701.433 0.590958 0.0293143 +701.756 0.599504 0.0344613 +702.08 0.601895 0.0322661 +702.403 0.612118 0.0311078 +702.726 0.593586 0.0308482 +703.05 0.59342 0.032416 +703.373 0.586261 0.0314058 +703.696 0.585284 0.0330124 +704.019 0.614425 0.0296919 +704.342 0.604172 0.0341668 +704.665 0.604788 0.0303507 +704.988 0.602451 0.0328664 +705.311 0.569533 0.0348221 +705.634 0.562828 0.029027 +705.957 0.572027 0.033083 +706.28 0.574146 0.0308819 +706.602 0.598662 0.0290074 +706.925 0.600156 0.0307929 +707.247 0.596965 0.0311694 +707.57 0.584376 0.0307132 +707.892 0.583197 0.028681 +708.215 0.594607 0.0294714 +708.537 0.570324 0.0311811 +708.859 0.564883 0.0287208 +709.182 0.570181 0.0312105 +709.504 0.566761 0.0306581 +709.826 0.589028 0.0312426 +710.148 0.589582 0.0304727 +710.47 0.592177 0.0298808 +710.792 0.583482 0.0296105 +711.114 0.587726 0.0299413 +711.436 0.586639 0.0307705 +711.758 0.574057 0.0313847 +712.079 0.591464 0.0324114 +712.401 0.577534 0.0299114 +712.723 0.583448 0.0282851 +713.044 0.576672 0.0325902 +713.366 0.589227 0.0287495 +713.687 0.594011 0.0302469 +714.009 0.589903 0.0291838 +714.33 0.589356 0.0292225 +714.652 0.583185 0.0303942 +714.973 0.593902 0.0285931 +715.294 0.58153 0.0307238 +715.615 0.586407 0.0316067 +715.936 0.59576 0.030383 +716.257 0.578766 0.0330603 +716.578 0.600498 0.0319397 +716.899 0.60386 0.0339516 +717.22 0.587555 0.0309088 +717.541 0.601989 0.0300648 +717.862 0.598223 0.0310932 +718.182 0.58935 0.0325868 +718.503 0.575619 0.0298679 +718.824 0.599702 0.0319173 +719.144 0.596405 0.0295792 +719.465 0.579453 0.0310685 +719.785 0.581857 0.0354539 +720.105 0.598952 0.0297556 +720.426 0.578401 0.0333904 +720.746 0.59302 0.0330676 +721.066 0.588245 0.0346719 +721.386 0.579071 0.0317136 +721.706 0.597385 0.0297824 +722.026 0.592311 0.0335772 +722.346 0.58317 0.0326709 +722.666 0.589741 0.036895 +722.986 0.587071 0.033553 +723.306 0.590615 0.0308245 +723.626 0.587175 0.0317767 +723.946 0.589365 0.0317444 +724.265 0.572202 0.0366396 +724.585 0.572484 0.0380243 +724.904 0.583332 0.0317413 +725.224 0.58596 0.0301587 +725.543 0.573233 0.0375397 +725.863 0.575455 0.0328122 +726.182 0.588139 0.0304016 +726.501 0.567124 0.0341404 +726.82 0.575276 0.0334414 +727.14 0.567261 0.0319789 +727.459 0.563214 0.0329479 +727.778 0.568435 0.0326287 +728.097 0.567714 0.0330952 +728.416 0.560075 0.0334835 +728.734 0.561562 0.0350045 +729.053 0.575493 0.0339702 +729.372 0.5637 0.0347561 +729.691 0.544671 0.030025 +730.009 0.563283 0.0346777 +730.328 0.555946 0.0344928 +730.647 0.571064 0.0331103 +730.965 0.570525 0.0330407 +731.283 0.565794 0.0299156 +731.602 0.54797 0.0295002 +731.92 0.538003 0.0332292 +732.238 0.540028 0.0328514 +732.557 0.560395 0.0338116 +732.875 0.542079 0.0356155 +733.193 0.532878 0.0327451 +733.511 0.558343 0.0354437 +733.829 0.570824 0.0299829 +734.147 0.574871 0.0333644 +734.465 0.548686 0.0303712 +734.782 0.53766 0.0343189 +735.1 0.542772 0.0319307 +735.418 0.546541 0.0334878 +735.735 0.555294 0.0305905 +736.053 0.543992 0.0360653 +736.371 0.552862 0.0326146 +736.688 0.556962 0.036478 +737.005 0.56719 0.032898 +737.323 0.545023 0.0333255 +737.64 0.541798 0.0314708 +737.957 0.531255 0.0400537 +738.274 0.531023 0.034583 +738.592 0.553722 0.034424 +738.909 0.576907 0.0336095 +739.226 0.564575 0.0313488 +739.543 0.551756 0.0309862 +739.86 0.568188 0.0322115 +740.176 0.570694 0.0321682 +740.493 0.555552 0.0424407 +740.81 0.542445 0.0376869 +741.127 0.553418 0.0332199 +741.443 0.529646 0.031511 +741.76 0.549652 0.0358277 +742.076 0.558962 0.0384082 +742.393 0.571116 0.0314491 +742.709 0.580253 0.0330698 +743.025 0.572058 0.0352556 +743.342 0.558659 0.0343983 +743.658 0.577667 0.0326625 +743.974 0.546902 0.0325538 +744.29 0.576221 0.0373525 +744.606 0.572196 0.038441 +744.922 0.551113 0.0334145 +745.238 0.555017 0.0315548 +745.554 0.564335 0.030366 +745.87 0.550381 0.0377499 +746.186 0.54471 0.0358074 +746.501 0.544372 0.0310092 +746.817 0.560815 0.0336819 +747.132 0.58055 0.0327317 +747.448 0.548818 0.0320504 +747.763 0.537438 0.0300465 +748.079 0.518969 0.0360061 +748.394 0.54048 0.0298564 +748.71 0.532809 0.0315933 +749.025 0.522945 0.0342235 +749.34 0.536727 0.0346545 +749.655 0.540802 0.0318589 +749.97 0.544994 0.0286103 +750.285 0.515331 0.0312239 +750.6 0.509121 0.0302846 +750.915 0.519378 0.0336343 +751.23 0.524415 0.0309258 +751.545 0.536222 0.02937 +751.859 0.522533 0.0311989 +752.174 0.509593 0.028469 +752.489 0.497657 0.0309337 +752.803 0.504869 0.0319697 +753.118 0.508398 0.0309006 +753.432 0.520047 0.0298089 +753.747 0.518769 0.031711 +754.061 0.499381 0.0324986 +754.375 0.494122 0.0315297 +754.689 0.500692 0.0310585 +755.004 0.500584 0.0299584 +755.318 0.50252 0.0341313 +755.632 0.48521 0.0343546 +755.946 0.494647 0.0370704 +756.26 0.47894 0.0316463 +756.573 0.482944 0.0303157 +756.887 0.478972 0.0293221 +757.201 0.48742 0.0333525 +757.515 0.495809 0.0291085 +757.828 0.481962 0.0288675 +758.142 0.466289 0.0300821 +758.455 0.463702 0.0332455 +758.769 0.472154 0.0311997 +759.082 0.462343 0.0333868 +759.396 0.464223 0.0305511 +759.709 0.475432 0.0302095 +760.022 0.48332 0.0301545 +760.335 0.480461 0.0289192 +760.648 0.47579 0.029441 +760.962 0.480171 0.0282019 +761.275 0.461625 0.028014 +761.587 0.46571 0.0335963 +761.9 0.476812 0.0309805 +762.213 0.471196 0.0311535 +762.526 0.477969 0.0288721 +762.839 0.47745 0.0306232 +763.151 0.473731 0.0299615 +763.464 0.459988 0.0323478 +763.776 0.458455 0.0302171 +764.089 0.469148 0.031061 +764.401 0.469795 0.0298841 +764.714 0.468504 0.0283195 +765.026 0.4664 0.0296671 +765.338 0.46865 0.0285082 +765.65 0.471636 0.0309584 +765.962 0.467753 0.0288404 +766.275 0.470096 0.0284387 +766.587 0.457026 0.0323591 +766.898 0.45588 0.0302618 +767.21 0.466277 0.0330532 +767.522 0.461446 0.0334344 +767.834 0.472927 0.0315779 +768.146 0.474462 0.0297973 +768.457 0.491338 0.0290766 +768.769 0.469495 0.0288502 +769.08 0.46175 0.0303557 +769.392 0.460135 0.0297195 +769.703 0.45859 0.0334191 +770.015 0.460849 0.0332581 +770.326 0.463733 0.0292005 +770.637 0.474827 0.0316299 +770.948 0.471648 0.0280418 +771.26 0.457737 0.0285002 +771.571 0.464498 0.0295879 +771.882 0.44985 0.0301881 +772.193 0.471399 0.0315716 +772.504 0.477654 0.030284 +772.814 0.473344 0.0293881 +773.125 0.475947 0.0307987 +773.436 0.479723 0.0287733 +773.747 0.484213 0.0315456 +774.057 0.495681 0.0292686 +774.368 0.497591 0.0335146 +774.678 0.506496 0.0342109 +774.989 0.493692 0.0336961 +775.299 0.512141 0.0332503 +775.609 0.493281 0.0315039 +775.919 0.50613 0.029808 +776.23 0.518301 0.0402589 +776.54 0.504315 0.0314382 +776.85 0.502739 0.034689 +777.16 0.496344 0.0341912 +777.47 0.498381 0.0317182 +777.78 0.499185 0.0312529 +778.089 0.487674 0.0318807 +778.399 0.490931 0.0284584 +778.709 0.48376 0.030613 +779.019 0.498135 0.0295682 +779.328 0.480746 0.0315728 +779.638 0.462388 0.0298491 +779.947 0.453142 0.0324409 +780.257 0.456875 0.0331757 +780.566 0.462936 0.029717 +780.875 0.455526 0.0298309 +781.185 0.463686 0.0297154 +781.494 0.465423 0.0300218 +781.803 0.4593 0.0289029 +782.112 0.458264 0.0289265 +782.421 0.468068 0.0298886 +782.73 0.460669 0.0328707 +783.039 0.481828 0.0298605 +783.348 0.483044 0.0311964 +783.656 0.486285 0.0324424 +783.965 0.468356 0.0302964 +784.274 0.479231 0.0340396 +784.582 0.480129 0.0326589 +784.891 0.476129 0.0299666 +785.199 0.473229 0.0301307 +785.508 0.475349 0.0305231 +785.816 0.508045 0.0299997 +786.124 0.475841 0.0321033 +786.432 0.461897 0.0298769 +786.741 0.482641 0.0324359 +787.049 0.474969 0.0305666 +787.357 0.471938 0.0301869 +787.665 0.464569 0.0300639 +787.973 0.470885 0.031123 +788.281 0.479502 0.031125 +788.588 0.474279 0.0305469 +788.896 0.477908 0.0329186 +789.204 0.483218 0.0329043 +789.511 0.481965 0.0338194 +789.819 0.49793 0.0336943 +790.126 0.495456 0.0320863 +790.434 0.517704 0.0317677 +790.741 0.496246 0.037687 +791.049 0.509227 0.0335077 +791.356 0.51495 0.0316888 +791.663 0.509734 0.0379308 +791.97 0.517525 0.0362853 +792.277 0.510078 0.034211 +792.584 0.483866 0.0379795 +792.891 0.467553 0.0300293 +793.198 0.469822 0.0331605 +793.505 0.471915 0.0313217 +793.812 0.457292 0.0300799 +794.119 0.447325 0.0292712 +794.425 0.460034 0.0313977 +794.732 0.457401 0.0318315 +795.038 0.502939 0.0304191 +795.345 0.492365 0.033166 +795.651 0.486079 0.0329532 +795.958 0.501696 0.0356552 +796.264 0.480498 0.0362322 +796.57 0.489355 0.0314284 +796.876 0.509117 0.0352206 +797.183 0.537177 0.0313034 +797.489 0.528045 0.0342307 +797.795 0.527707 0.0315038 +798.101 0.538203 0.034152 +798.406 0.554143 0.037344 +798.712 0.542333 0.0323327 +799.018 0.550233 0.036784 +799.324 0.524045 0.0320027 +799.629 0.548619 0.0341161 +799.935 0.530901 0.0335224 +800.24 0.576676 0.0339562 +800.546 0.582409 0.0410324 +800.851 0.544801 0.0359171 +801.157 0.561008 0.0384545 +801.462 0.580871 0.0362456 +801.767 0.572528 0.0359954 +802.072 0.57711 0.0347733 +802.377 0.567713 0.0367473 +802.682 0.551787 0.0321518 +802.987 0.557532 0.0374502 +803.292 0.553483 0.0380711 +803.597 0.563631 0.0357474 +803.902 0.556487 0.0376004 +804.207 0.564993 0.0421552 +804.511 0.5541 0.0355413 +804.816 0.539148 0.0374116 +805.12 0.502262 0.032716 +805.425 0.497826 0.0337669 +805.729 0.519177 0.0325661 +806.034 0.497712 0.0302748 +806.338 0.522508 0.0335725 +806.642 0.503648 0.032304 +806.946 0.51072 0.0384972 +807.251 0.540532 0.0319409 +807.555 0.491738 0.0314612 +807.859 0.539879 0.0323729 +808.162 0.542829 0.034929 +808.466 0.545231 0.0379279 +808.77 0.537416 0.0344004 +809.074 0.53789 0.0365726 +809.378 0.571448 0.0359715 +809.681 0.553084 0.0479653 +809.985 0.568987 0.035978 +810.288 0.611285 0.0421301 +810.592 0.558276 0.0346725 +810.895 0.569606 0.0419553 +811.198 0.592181 0.042674 +811.502 0.605605 0.0437411 +811.805 0.636961 0.0386894 +812.108 0.580279 0.0424028 +812.411 0.561413 0.0367032 +812.714 0.567239 0.0376454 +813.017 0.547089 0.0328613 +813.32 0.55294 0.0322086 +813.623 0.553474 0.0348743 +813.926 0.585483 0.0363626 +814.228 0.59519 0.0360118 +814.531 0.611724 0.0386641 +814.833 0.606006 0.0357493 +815.136 0.600037 0.0391627 +815.438 0.605367 0.0398666 +815.741 0.586403 0.0383424 +816.043 0.574048 0.0355256 +816.345 0.586923 0.035907 +816.648 0.590784 0.0395323 +816.95 0.607444 0.0404561 +817.252 0.623491 0.0408207 +817.554 0.597015 0.0461405 +817.856 0.622788 0.0399739 +818.158 0.624035 0.0390888 +818.46 0.644078 0.0409552 +818.761 0.658839 0.0428061 +819.063 0.612247 0.0381264 +819.365 0.604548 0.0384584 +819.666 0.583159 0.0350323 +819.968 0.582579 0.0369943 +820.269 0.56511 0.0384575 +820.571 0.599095 0.0410713 +820.872 0.59083 0.0442579 +821.173 0.596848 0.0379478 +821.475 0.623419 0.0368687 +821.776 0.615436 0.0497453 +822.077 0.609811 0.041157 +822.378 0.647384 0.0415349 +822.679 0.657904 0.0466911 +822.98 0.629769 0.0446454 +823.281 0.646234 0.0413733 +823.581 0.644712 0.0393263 +823.882 0.658973 0.0406557 +824.183 0.644684 0.0420204 +824.483 0.669977 0.0413001 +824.784 0.588804 0.0460724 +825.084 0.627286 0.03698 +825.385 0.639062 0.0423538 +825.685 0.649686 0.047504 +825.985 0.615921 0.0415616 +826.286 0.618212 0.0364799 +826.586 0.599358 0.0398824 +826.886 0.625471 0.044072 +827.186 0.607821 0.0396996 +827.486 0.659164 0.0382515 +827.786 0.639591 0.0402897 +828.086 0.670614 0.0457159 +828.386 0.626103 0.0366217 +828.685 0.622889 0.0399933 +828.985 0.645313 0.0406605 +829.285 0.665852 0.0413793 +829.584 0.667006 0.0405444 +829.884 0.68221 0.0408422 +830.183 0.701672 0.0420037 +830.482 0.663506 0.0439649 +830.782 0.687668 0.0411885 +831.081 0.726174 0.0422757 +831.38 0.710294 0.0418266 +831.679 0.723683 0.04232 +831.978 0.725205 0.0469256 +832.277 0.699759 0.0404471 +832.576 0.661313 0.0400649 +832.875 0.678147 0.0395181 +833.174 0.678373 0.0408175 +833.472 0.669932 0.0457291 +833.771 0.72679 0.0500605 +834.07 0.645925 0.0398786 +834.368 0.702045 0.0498931 +834.667 0.726601 0.0464165 +834.965 0.734784 0.0518451 +835.263 0.727552 0.0443727 +835.562 0.693597 0.0493369 +835.86 0.687644 0.0405149 +836.158 0.671701 0.0397927 +836.456 0.695346 0.044312 +836.754 0.743261 0.0443395 +837.052 0.76107 0.0562341 +837.35 0.756271 0.049257 +837.648 0.7518 0.0462969 +837.946 0.749951 0.0462109 +838.243 0.7729 0.0536236 +838.541 0.753742 0.0528466 +838.839 0.734552 0.0443817 +839.136 0.761371 0.0458231 +839.434 0.729477 0.0410123 +839.731 0.73896 0.051427 +840.028 0.760926 0.0449277 +840.326 0.690731 0.0539243 +840.623 0.727856 0.049651 +840.92 0.759805 0.0470471 +841.217 0.76575 0.050162 +841.514 0.789345 0.0455834 +841.811 0.725763 0.0503215 +842.108 0.825354 0.0449185 +842.405 0.807558 0.0489371 +842.701 0.77625 0.0502582 +842.998 0.744803 0.047018 +843.295 0.803966 0.0465682 +843.591 0.762032 0.0526688 +843.888 0.768946 0.0583749 +844.184 0.748223 0.0453922 +844.481 0.738664 0.045964 +844.777 0.777956 0.0448096 +845.073 0.751849 0.0465908 +845.369 0.775199 0.0468082 +845.665 0.82189 0.0503622 +845.962 0.778715 0.0464545 +846.258 0.81243 0.0507204 +846.553 0.814326 0.0501506 +846.849 0.75504 0.0445088 +847.145 0.712489 0.0489688 +847.441 0.770972 0.0630369 +847.737 0.797241 0.0469129 +848.032 0.744711 0.0435029 +848.328 0.776843 0.0488394 +848.623 0.760924 0.0542879 +848.919 0.767533 0.0457016 +849.214 0.771708 0.052886 +849.509 0.786281 0.0466834 +849.805 0.818148 0.0525345 +850.1 0.843293 0.051857 +850.395 0.864676 0.0666683 +850.69 0.843946 0.0501806 +850.985 0.790672 0.0593455 +851.28 0.780452 0.0436372 +851.575 0.786924 0.054512 +851.869 0.823576 0.0486607 +852.164 0.743056 0.0560804 +852.459 0.756835 0.0452691 +852.753 0.779636 0.0499919 +853.048 0.666395 0.0436665 +853.342 0.68044 0.0466099 +853.637 0.703585 0.0457622 +853.931 0.696424 0.041391 +854.225 0.665451 0.0437807 +854.52 0.74711 0.0497545 +854.814 0.715503 0.0493889 +855.108 0.693536 0.0401881 +855.402 0.717595 0.0485217 +855.696 0.734491 0.0449131 +855.99 0.743449 0.0458442 +856.283 0.745841 0.0418643 +856.577 0.779059 0.0465935 +856.871 0.745798 0.0472604 +857.164 0.756101 0.0491214 +857.458 0.827325 0.045297 +857.751 0.811695 0.0481812 +858.045 0.735561 0.0586892 +858.338 0.772291 0.0434434 +858.632 0.714342 0.0465045 +858.925 0.758145 0.0485414 +859.218 0.756607 0.0459267 +859.511 0.835746 0.0489136 +859.804 0.794034 0.0509284 +860.097 0.763727 0.0539569 +860.39 0.797337 0.0474293 +860.683 0.83393 0.0493979 +860.976 0.808106 0.0463319 +861.268 0.830624 0.04619 +861.561 0.818045 0.0532449 +861.854 0.950802 0.0524534 +862.146 0.826856 0.0537291 +862.439 0.798768 0.0455337 +862.731 0.753283 0.0587547 +863.023 0.808218 0.0473111 +863.316 0.8079 0.0631359 +863.608 0.823477 0.0518033 +863.9 0.759566 0.0450218 +864.192 0.756049 0.044786 +864.484 0.813976 0.0666071 +864.776 0.809196 0.0518238 +865.068 0.807301 0.0485325 +865.36 0.783184 0.0466981 +865.651 0.853318 0.0587228 +865.943 0.773988 0.0454597 +866.235 0.823625 0.052128 +866.526 0.836859 0.0506749 +866.818 0.800318 0.0457736 +867.109 0.815579 0.0532148 +867.4 0.820273 0.0456343 +867.692 0.796014 0.0548978 +867.983 0.758279 0.0566296 +868.274 0.840299 0.0462718 +868.565 0.847685 0.0545165 +868.856 0.833204 0.0529598 +869.147 0.822322 0.0474649 +869.438 0.834472 0.0489846 +869.729 0.835121 0.0584197 +870.02 0.848958 0.0646508 +870.31 0.86086 0.0552112 +870.601 0.797363 0.0568874 +870.892 0.814385 0.0524096 +871.182 0.775926 0.0481121 +871.473 0.854375 0.0467963 +871.763 0.816242 0.0472453 +872.053 0.850469 0.0475086 +872.343 0.807151 0.0562472 +872.634 0.793155 0.0490103 +872.924 0.799976 0.0506747 +873.214 0.792468 0.0486992 +873.504 0.811375 0.0464993 +873.794 0.86456 0.06991 +874.084 0.872162 0.0501665 +874.373 0.825339 0.050673 +874.663 0.789072 0.0503136 +874.953 0.83129 0.0497175 +875.242 0.848539 0.0536462 +875.532 0.803103 0.0506352 +875.821 0.805738 0.0494973 +876.111 0.792315 0.0491531 +876.4 0.821777 0.0544084 +876.689 0.82569 0.053746 +876.978 0.807034 0.0549065 +877.268 0.819824 0.0514371 +877.557 0.839354 0.0505359 +877.846 0.868996 0.0489987 +878.134 0.832026 0.055442 +878.423 0.832179 0.055447 +878.712 0.832332 0.055452 diff --git a/include/qsimDetectorConstruction.hh b/include/qsimDetectorConstruction.hh index 1ce60f3..caaf704 100644 --- a/include/qsimDetectorConstruction.hh +++ b/include/qsimDetectorConstruction.hh @@ -1,10 +1,12 @@ - +//From brads_qsim branch #ifndef qsimDetectorConstruction_h #define qsimDetectorConstruction_h 1 #include "globals.hh" +#include "G4SubtractionSolid.hh" #include "G4VUserDetectorConstruction.hh" - +#include +#include //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... class qsimDetectorConstruction : public G4VUserDetectorConstruction @@ -12,19 +14,22 @@ class qsimDetectorConstruction : public G4VUserDetectorConstruction public: qsimDetectorConstruction(); ~qsimDetectorConstruction(); - //void StandModeSet(); - void DetModeSet(G4int ); - void StandModeSet(G4int ); - public: + void DetModeSet(G4int ); + void StandModeSet(G4int ); + void LGReflectivitySet(G4double ); + public: G4VPhysicalVolume* Construct(); private: + G4int fDetMode; + G4int fStandMode; + G4double det_x; + G4double det_y; + G4double det_z; + G4double quartz_x; G4double quartz_y; G4double quartz_z; - //G4int fStandMode; - G4int fDetMode; - G4int fStandMode; G4double quartz_zPos; @@ -40,12 +45,17 @@ class qsimDetectorConstruction : public G4VUserDetectorConstruction G4double rout; G4double lngth; - public: - G4double fDetAngle, fQuartzPolish; - // POSSCAN - G4double fDetPosX, fDetPosY; - - + std::ifstream myfile; + std::ofstream outfile; + std::ifstream myfile2; + + public: + G4int gasType; + G4double beam_angle; + G4double fDetAngle, fQuartzPolish; + G4double fLGReflectivity; + // POSSCAN + G4double fDetPosX, fDetPosY; }; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/include/qsimMessenger.hh b/include/qsimMessenger.hh index 84a2aab..bd7b875 100644 --- a/include/qsimMessenger.hh +++ b/include/qsimMessenger.hh @@ -16,6 +16,7 @@ class qsimDetectorConstruction; class qsimEventAction; class qsimPrimaryGeneratorAction; class qsimSteppingAction; +class qsimOpticalPhysics; class G4UIcmdWithAnInteger; class G4UIcmdWithADoubleAndUnit; @@ -31,6 +32,7 @@ class qsimMessenger : public G4UImessenger { void SetIO( qsimIO *io ){ fIO = io; } void SetPriGen( qsimPrimaryGeneratorAction *pg ){ fprigen = pg; } + void SetOptPhys( qsimOpticalPhysics *op ){ foptical = op; } void SetDetCon( qsimDetectorConstruction *dc ){ fdetcon= dc; } void SetEvAct( qsimEventAction *ev ){ fevact = ev; } void SetStepAct( qsimSteppingAction *st ){ fStepAct = st; } @@ -43,16 +45,20 @@ class qsimMessenger : public G4UImessenger { qsimEventAction *fevact; qsimPrimaryGeneratorAction *fprigen; qsimSteppingAction *fStepAct; + qsimOpticalPhysics *foptical; G4UIdirectory *fRemollDir; G4UIcmdWithAnInteger *seedCmd; G4UIcmdWithAString *fileCmd; // + G4UIcmdWithABool *fAllowCerenkovCmd; + G4UIcmdWithABool *fAllowScintillationCmd; //G4UIcmdWithAnInteger *fStandModeCmd; G4UIcmdWithAnInteger *fDetModeCmd; G4UIcmdWithAnInteger *fStandModeCmd; G4UIcmdWithAnInteger *fSourceModeCmd; + G4UIcmdWithADouble *fLGReflectivityCmd; G4UIcmdWithADouble *fQuartzPolishCmd; G4UIcmdWithADoubleAndUnit *fDetAngleCmd; // POSSCAN diff --git a/include/qsimOpticalPhysics.hh b/include/qsimOpticalPhysics.hh index 6dfe4fb..6db17fb 100644 --- a/include/qsimOpticalPhysics.hh +++ b/include/qsimOpticalPhysics.hh @@ -18,8 +18,10 @@ class qsimOpticalPhysics : public G4VPhysicsConstructor { public: + qsimOpticalPhysics(); qsimOpticalPhysics(G4bool toggle=true); - virtual ~qsimOpticalPhysics(); + ~qsimOpticalPhysics(); + //virtual ~qsimOpticalPhysics(); virtual void ConstructParticle(); virtual void ConstructProcess(); @@ -32,6 +34,9 @@ class qsimOpticalPhysics : public G4VPhysicsConstructor G4OpMieHG* GetMieHGScatteringProcess() {return theMieHGScatteringProcess;} G4OpBoundaryProcess* GetBoundaryProcess() { return theBoundaryProcess;} + void AllowCerenkovSet(G4bool ); + void AllowScintillationSet(G4bool ); + void SetNbOfPhotonsCerenkov(G4int); private: @@ -43,6 +48,9 @@ private: G4OpRayleigh* theRayleighScattering; G4OpMieHG* theMieHGScatteringProcess; G4OpBoundaryProcess* theBoundaryProcess; + + bool fAllowCerenkov; + bool fAllowScintillation; G4bool AbsorptionOn; diff --git a/macros/get_pe.C b/macros/get_pe.C new file mode 100644 index 0000000..ccdc329 --- /dev/null +++ b/macros/get_pe.C @@ -0,0 +1,351 @@ +#include +#include + +using namespace std; + + /* +const double w = 0.5473; +const double Q0 = 73.91; +const double sigma_0 = 0.5788; +const double alpha = 0.06981; +const double mu = 5.871; +const double Q1 = 25.4; +const double sigma_1 = 12.15; + +//Sub functions +static const double degtorad = TMath::Pi()/180.; + +static const double twopi = 2*TMath::Pi(); + +Double_t fit_func_yxzhao_pe(Double_t *x, Double_t *par); +Double_t fit_func_yxzhao_bg(Double_t *x, Double_t *par); +*/ + +int det_id=2; //1 for the upstream, 2 for the downstream + +void get_pe(std::string fname = "qsim.root", double user_reflectivity = 0.9, double user_cerenkov = 1.0, double user_scintillation = 1.0, double user_z_pos = -11.0){ + + TFile *root_file = new TFile(fname.c_str()); + TTree *T = (TTree*)root_file->Get("T"); + + int hit_n=0; + int hit_det[100000]={0}; + int hit_pid[100000]={0}; + double hit_x[100000] = {0}; + double hit_z[100000] = {0}; + double hit_e[100000] = {0}; + + double hit_y[100000] = {0}; + + + T->SetBranchAddress("hit.n",&hit_n); + T->SetBranchAddress("hit.det",hit_det); + T->SetBranchAddress("hit.pid",hit_pid); + T->SetBranchAddress("hit.x", hit_x); + T->SetBranchAddress("hit.y", hit_y); + T->SetBranchAddress("hit.z", hit_z); + T->SetBranchAddress("hit.e", hit_e); + + long N_entries=T->GetEntries(); + + TH1F *h_npe=new TH1F("npe","number of p.e",30,0,30); + TH1F *h_npeN=new TH1F("npeN","Normalized number of p.e",30,0,30); + TH1F *h = new TH1F("h", "Photon Energy", 10000, 0, 10); + //#TH2F *h_hp = new TH2F("hp", "Hit position", 1000, -.1, .1, 1000, .2, .4); + //TH1F *h_QDC_sim=new TH1F("h_QDC_sim","smeared G4 simulation in QDC channel",4096,0,4096); + + int bins[100000] = {0}; + int sum = 0; + + for(long i=0; iGetEntry(i); + //cout<Fill(hit_x[j], hit_z[j]); + npe++; + sum++; + //h->Fill(hit_e[j] / (10^-9)); + // } + } + } +/* + // This is the smear + double parameters[8]={w,Q0, sigma_0,alpha,mu, Q1, sigma_1, npe }; + + TF1 *func_gen; + + if(npe==0){ + func_gen=new TF1("func_gen",fit_func_yxzhao_bg,0,4096,5); + func_gen->SetParameters(parameters); + }else{ + func_gen=new TF1("func_gen", fit_func_yxzhao_pe, 0, 4096, 8); + func_gen->SetParameters(parameters); + } + + double QDC_chan_gen=func_gen->GetRandom(); + cout<<"##################### "<Fill(QDC_chan_gen); + + delete func_gen; +*/ + bins[i] = npe; + } + + cout<Fill(bins[k]); + h_npeN->Fill(bins[k], 1.0/sum); + } + //h_QDC_sim->SetAxisRange(68, 600); + + TCanvas * c1 = new TCanvas(); + c1->Divide(2,1); + c1->cd(1); + h_npe->Draw(); + gPad->SetLogy(); + c1->cd(2); + h_npeN->Draw(); + gPad->SetLogy(); + std::cout << "Mean PE/event = " << h_npe->GetMean() << ", RMS = " << h_npe->GetRMS() << ", and resolution = " << h_npe->GetRMS()/h_npe->GetMean() << " for " << N_entries << " entries" << std::endl; + + std::string fnameOutPDF = fname.substr(0,fname.find(".root")) + ".pdf"; + c1->SaveAs(fnameOutPDF.c_str()); + + // Save data + double ref_x_pos = atof((fname.substr(fname.find("ees_")+4,fname.find("_x")-fname.find("ees_")-4)).c_str()); + double ref_angle = atof((fname.substr(fname.find("out_")+4,fname.find("_degrees")-fname.find("out_")-4)).c_str()); + double ref_reflectivity = user_reflectivity; + double ref_cerenkov = user_cerenkov; + double ref_scintillation = user_scintillation; + double ref_z_pos = user_z_pos; + + double oldx_pos = 0.0; + double oldangle = 0.0; + double oldavg = 0.0; + double oldavg_err = 0.0; + double oldrms = 0.0; + double oldrms_err = 0.0; + double oldres = 0.0; + double oldN_en = 0.0; + double oldreflectivity = 0.0; + double oldcerenkov = 0.0; + double oldscintillation = 0.0; + double oldz_pos = 0.0; + + double x_pos = 0.0; + double angle = 0.0; + double avg = 0.0; + double avg_err = 0.0; + double rms = 0.0; + double rms_err = 0.0; + double res = 0.0; + double N_en = 0.0; + double reflectivity = 0.0; + double cerenkov = 0.0; + double scintillation = 0.0; + double z_pos = 0.0; + + std::cout << "X = " << ref_x_pos << ", angle = " << ref_angle <AccessPathName("scans.root")) { + // Old file exists, read it and add new entries + old_file = TFile::Open("scans.root"); + old_file->GetObject("scans", oldtree); + new_file.cd(); + if (!oldtree) { + std::cout << "ERROR: Dead scans tree" ; + return; + } + newtree = oldtree->CloneTree(0); + int nent = oldtree->GetEntries(); + +// TLeaf* angleL = oldtree->GetLeaf("angle"); +// TLeaf* x_posL = oldtree->GetLeaf("x_pos"); + + // Clear out prior instance if exists + bool prior = true; + if (oldtree->GetBranch("reflectivity")) { + oldtree->SetBranchAddress("reflectivity",&oldreflectivity); + newtree->SetBranchAddress("reflectivity",&reflectivity); + } + else { + newtree->Branch("reflectivity",&reflectivity); + prior = false; + } + if (oldtree->GetBranch("cerenkov")) { + oldtree->SetBranchAddress("cerenkov",&oldcerenkov); + newtree->SetBranchAddress("cerenkov",&cerenkov); + } + else { + newtree->Branch("cerenkov",&cerenkov); + prior = false; + } + if (oldtree->GetBranch("scintillation")) { + oldtree->SetBranchAddress("scintillation",&oldscintillation); + newtree->SetBranchAddress("scintillation",&scintillation); + } + else { + newtree->Branch("scintillation",&scintillation); + prior = false; + } + if (oldtree->GetBranch("z_pos")) { + oldtree->SetBranchAddress("z_pos",&oldz_pos); + newtree->SetBranchAddress("z_pos",&z_pos); + } + else { + newtree->Branch("z_pos",&z_pos); + prior = false; + } + oldtree->SetBranchAddress("angle",&oldangle); + oldtree->SetBranchAddress("x_pos",&oldx_pos); + oldtree->SetBranchAddress("avg_pes",&oldavg); + oldtree->SetBranchAddress("avg_pes_err",&oldavg_err); + oldtree->SetBranchAddress("rms_pes",&oldrms); + oldtree->SetBranchAddress("rms_pes_err",&oldrms_err); + oldtree->SetBranchAddress("res",&oldres); + oldtree->SetBranchAddress("nentries",&oldN_en); + newtree->SetBranchAddress("angle",&angle); + newtree->SetBranchAddress("x_pos",&x_pos); + newtree->SetBranchAddress("avg_pes",&avg); + newtree->SetBranchAddress("avg_pes_err",&avg_err); + newtree->SetBranchAddress("rms_pes",&rms); + newtree->SetBranchAddress("rms_pes_err",&rms_err); + newtree->SetBranchAddress("res",&res); + newtree->SetBranchAddress("nentries",&N_en); + for (auto j : ROOT::TSeqI(nent) ) { +// x_posL->GetBranch()->GetEntry(j); +// angleL->GetBranch()->GetEntry(j); + oldtree->GetEntry(j); + + if (ref_x_pos == oldx_pos && ref_angle == oldangle && (!prior || (ref_reflectivity == oldreflectivity && ref_cerenkov == oldcerenkov && ref_scintillation == oldscintillation && ref_z_pos == oldz_pos))) { + //if (ref_x_pos == x_posL->GetValue() && ref_angle == angleL->GetValue()) + std::cout << "TEST 1" << std::endl; + continue; + } + x_pos = oldx_pos; + angle = oldangle; + avg = oldavg; + avg_err = oldavg_err; + rms = oldrms; + rms_err = oldrms_err; + res = oldres; + N_en = oldN_en; + reflectivity = oldreflectivity; + cerenkov = oldcerenkov; + scintillation = oldscintillation; + z_pos = oldz_pos; + if (!oldtree->GetBranch("reflectivity")) { + reflectivity = 0.9; + } + if (!oldtree->GetBranch("cerenkov")) { + cerenkov = 1.0; + } + if (!oldtree->GetBranch("scintillation")) { + scintillation = 1.0; + } + if (!oldtree->GetBranch("z_pos")) { + z_pos = -11.0; + } + newtree->Fill(); + } + + // Append current run to end + old_file->Close(); + gSystem->Exec("rm scans.root"); + delete old_file; + } + else { + // Old file doesn't exist, make a new one + new_file.cd(); + newtree = new TTree("scans","scans"); + + // Write new tree + newtree->Branch("angle",&angle); + newtree->Branch("x_pos",&x_pos); + newtree->Branch("avg_pes",&avg); + newtree->Branch("avg_pes_err",&avg_err); + newtree->Branch("rms_pes",&rms); + newtree->Branch("rms_pes_err",&rms_err); + newtree->Branch("res",&res); + newtree->Branch("nentries",&N_en); + newtree->Branch("reflectivity",&reflectivity); + newtree->Branch("cerenkov",&cerenkov); + newtree->Branch("scintillation",&scintillation); + newtree->Branch("z_pos",&z_pos); + } + + newtree->SetBranchAddress("angle",&angle); + newtree->SetBranchAddress("x_pos",&x_pos); + newtree->SetBranchAddress("avg_pes",&avg); + newtree->SetBranchAddress("avg_pes_err",&avg_err); + newtree->SetBranchAddress("rms_pes",&rms); + newtree->SetBranchAddress("rms_pes_err",&rms_err); + newtree->SetBranchAddress("res",&res); + newtree->SetBranchAddress("nentries",&N_en); + newtree->SetBranchAddress("reflectivity",&reflectivity); + newtree->SetBranchAddress("cerenkov",&cerenkov); + newtree->SetBranchAddress("scintillation",&scintillation); + newtree->SetBranchAddress("z_pos",&z_pos); + + x_pos = atof((fname.substr(fname.find("ees_")+4,fname.find("_x")-fname.find("ees_")-4)).c_str()); + angle = atof((fname.substr(fname.find("out_")+4,fname.find("_degrees")-fname.find("out_")-4)).c_str()); + avg = h_npe->GetMean(); + avg_err = h_npe->GetMeanError(); + rms = h_npe->GetRMS(); + rms_err = h_npe->GetRMSError(); + res = h_npe->GetRMS()/h_npe->GetMean(); + N_en = (double)N_entries; + + reflectivity = user_reflectivity; + cerenkov = user_cerenkov; + scintillation = user_scintillation; + z_pos = user_z_pos; + + std::cout << "TEST 2 X = " << x_pos << ", angle = " << angle <Fill(); + newtree->Write("scans",TObject::kOverwrite); + new_file.Close(); + + gSystem->Exec("mv localTmp.root scans.root"); + + //#h_hp->Draw(); +} +/* +//Fitting functions + +Double_t fit_func_yxzhao_pe(Double_t *x, Double_t *par){ + double xx = x[0]; + double s_real_sum=0; + double term_1 = xx-par[1]-par[0]/par[3]-par[7]*par[5]; + s_real_sum=TMath::PoissonI(par[7],par[4]) *( 1.0/par[6]/sqrt(twopi)/sqrt(par[7]) * exp(-1.0*pow + (term_1, 2)/2.0/par[7]/pow(par[6],2))); + return s_real_sum; +} + +Double_t fit_func_yxzhao_bg(Double_t *x, Double_t *par){ + double xx = x[0]; + double theta_f=0; + if(xx>=par[1]){ + theta_f=1.0; + }else{ + theta_f=0; + } + double background=exp(-1.0*par[4]) * ( (1-par[0])/par[2]/sqrt(twopi)*exp(-1.0*pow(xx-par[1],2)/2.0/pow(par[2],2) ) + par[0]*theta_f*par[3]*exp(-1.0*par[3]*(xx-par[1])) ); + + return background; +} */ diff --git a/macros/preserve.mac b/macros/preserve.mac new file mode 100644 index 0000000..1550038 --- /dev/null +++ b/macros/preserve.mac @@ -0,0 +1,31 @@ +#/vis/open OGLIQt 1200x600 +/qsim/fSourceMode 1 +/gun/particle e- + +/qsim/xmin 0.0 cm +/qsim/xmax 0.0 cm +/qsim/ymin -0.5 cm +/qsim/ymax 0.5 cm + +/qsim/z -11.0 cm + +/qsim/thetamin 0.0 deg +/qsim/thetamax 0.0 deg + +/qsim/fLGReflectivity 0.9 +/qsim/fAllowCerenkov 1 +/qsim/fAllowScintillation 1 + +#/gun/position 0.0 0.0 -11.0 cm +#/gun/direction 0 0 1 + +/qsim/emin 850 MeV +/qsim/emax 850 MeV + + +#/control/execute vis/vis.mac + +/run/initialize + +/qsim/filename qsimout_0.0_degrees_0.0_x.root +/run/beamOn 100000 diff --git a/macros/qsim.sh b/macros/qsim.sh new file mode 100755 index 0000000..7979080 --- /dev/null +++ b/macros/qsim.sh @@ -0,0 +1,114 @@ +#!/bin/bash + +# Takes 5 arguments +# +# Fixed value of non-scanned variable - default = 0.0 +# Variable to scan ("angle" or "x") - default = "angle" +# Min of scanned variable - default = -30 +# Max of scanned variable - default = 30 +# Step size - default = 0.5 + +if [ "$#" -lt 1 ] ; then + echo " ERROR, requires at least one input + " + echo " usage: ./qsim.sh fixed-non-scanned \"variable-to-scan\" min-of-scan (-30) max-of-scan (30) step-size (0.5) + Takes 9 arguments + + Fixed value of non-scanned variable - default = 0.0 + Variable to scan (\"angle\" or \"x\") - default = \"angle\" + Min of scanned variable - default = -30 + Max of scanned variable - default = 30 + Step size - default = 0.5 + Reflectivity of light guide - default = 0.9 + Cerenkov - default = 1.0 + Scintillation - default = 1.0 + z Position of beam origin - default = -11.0" + exit +fi +fixed="$1" +scanned="angle" +scanMin=-30.0 +scanMax=30.0 +scanStep=0.5 +reflectivity=0.9 +cerenkov=1 +scintillation=1 +#z_pos=-52.0 +z_pos=-11.0 +if [ "$#" -gt 1 ] ; then + scanned="$2" +fi +if [ "$#" -gt 2 ] ; then + scanMin=$3 +fi +if [ "$#" -gt 3 ] ; then + scanMax=$4 +fi +if [ "$#" -gt 4 ] ; then + scanStep=$5 +fi +if [ "$#" -gt 5 ] ; then + reflectivity=$6 +fi +if [ "$#" -gt 6 ] ; then + cerenkov=$(printf "%.0f" "$7") +fi +if [ "$#" -gt 7 ] ; then + scintillation=$(printf "%.0f" "$8") +fi +if [ "$#" -gt 8 ] ; then + z_pos=$9 +fi + +numSteps=$(printf "%.0f" "$(bc -l <<< \(${scanMax}-\(${scanMin}\)\)/$scanStep)") +fixed=$(printf "%.1f" "$(bc -l <<< 1.0*$fixed)") +z_point=-11.0 + +defaultName="qsimout.root" +outputName="qsimout.root" +name="0.0_degrees_0.0_x" + +for currentStep in `seq 0 $numSteps`; +do + point=$(printf "%.1f" "$(bc -l <<< \($scanMin+\(1.0*$currentStep*$scanStep\)\))") + if [[ "$scanned" == "angle" ]] ; then + name="${point}_degrees_${fixed}_x_${z_pos}_z_${reflectivity}_ref_${cerenkov}_cer_${scintillation}_scint" + cp preserve.mac sim_${name}.mac + sed -i 's;'"thetamin 0.0 deg"';'"thetamin ${point} deg"';g' sim_${name}.mac + sed -i 's;'"thetamax 0.0 deg"';'"thetamax ${point} deg"';g' sim_${name}.mac + sed -i 's;'"xmin 0.0 cm"';'"xmin ${fixed} cm"';g' sim_${name}.mac + sed -i 's;'"xmax 0.0 cm"';'"xmax ${fixed} cm"';g' sim_${name}.mac + # 0.1994 = sin(11.5 degrees), $fixed is the distance from 0.0, + is farther towards PMT, further back away from +z axis + z_point=$(printf "%.1f" "$(bc -l <<< ${z_pos}-\($fixed*0.1994\))") + #z_point=$(printf "%.1f" "$(bc -l <<< -11.0-\($fixed*0.1994\))") + fi + if [[ "$scanned" == "x" ]] ; then + name="${fixed}_degrees_${point}_x_${z_point}_z_${reflectivity}_ref_${cerenkov}_cer_${scintillation}_scint" + cp sim_0.0_degrees_0.0_x.mac sim_${name}.mac + sed -i 's;'"thetamin 0.0 deg"';'"thetamin ${fixed} deg"';g' sim_${name}.mac + sed -i 's;'"thetamax 0.0 deg"';'"thetamax ${fixed} deg"';g' sim_${name}.mac + sed -i 's;'"xmin 0.0 cm"';'"xmin ${point} cm"';g' sim_${name}.mac + sed -i 's;'"xmax 0.0 cm"';'"xmax ${point} cm"';g' sim_${name}.mac + # 0.1994 = sin(11.5 degrees), $point is the distance from 0.0, + is farther towards PMT, further back away from +z axis + z_point=$(printf "%.1f" "$(bc -l <<< ${z_pos}-\($point*0.1994\))") + #z_point=$(printf "%.1f" "$(bc -l <<< -11.0-\($point*0.1994\))") + fi + sed -i 's;'"fLGReflectivity 0.9"';'"fLGReflectivity ${reflectivity}"';g' sim_${name}.mac + sed -i 's;'"fAllowCerenkov 1"';'"fAllowCerenkov ${cerenkov}"';g' sim_${name}.mac + sed -i 's;'"fAllowScintillation 1"';'"fAllowScintillation ${scintillation}"';g' sim_${name}.mac + sed -i 's;'"z -11.0 cm"';'"z ${z_point} cm"';g' sim_${name}.mac + sed -i 's;'"0.0_degrees_0.0_x.root"';'"${name}.root"';g' sim_${name}.mac + cd ../ + source /share/apps/root-5.34.36-build/bin/thisroot.sh +# if [[ "$scintillation" != "1.0" ]] ; then +# ./build/qsim-cerenkovOnly macros/sim_${name}.mac +# elif [[ "$cerenkov" != "1.0" ]] ; then +# ./build/qsim-scintOnly macros/sim_${name}.mac +# else + ./build/qsim macros/sim_${name}.mac +# fi + cd - + source /share/apps/root-6.14.06-build/bin/thisroot.sh + root -l -b -q get_pe.C'("'../qsimout_${name}.root'",'${reflectivity}','${cerenkov}','${scintillation}','${z_point}')' + source /share/apps/root-5.34.36-build/bin/thisroot.sh +done diff --git a/macros/vis_brad.mac b/macros/vis_brad.mac new file mode 100644 index 0000000..4ef7f8f --- /dev/null +++ b/macros/vis_brad.mac @@ -0,0 +1,23 @@ +/vis/open OGLIQt 1200x600 + +/qsim/fSourceMode 1 +/gun/particle e- + +/qsim/xmin 0.0 cm +/qsim/xmax 0.0 cm +/qsim/ymin 0.0 cm +/qsim/ymax 0.0 cm + +/qsim/z -11.0 cm + +/qsim/thetamin 0.0 deg +/qsim/thetamax 0.0 deg + +#/gun/position 0.0 0.0 -11.0 cm +#/gun/direction 0 0 1 + +/qsim/emin 850 MeV +/qsim/emax 850 MeV + +/control/execute vis/vis.mac + diff --git a/qe.txt b/qe.txt new file mode 100644 index 0000000..d36223a --- /dev/null +++ b/qe.txt @@ -0,0 +1,548 @@ +160.01 0.672067 +160.02 0.970745 +160.03 1.26942 +160.04 1.56827 +160.049 1.86695 +163.338 2.16562 +163.347 2.4643 +163.357 2.76298 +163.367 3.06166 +163.376 3.36017 +163.386 3.65918 +166.674 3.95786 +166.684 4.25637 +166.694 4.55522 +166.704 4.85389 +166.713 5.15257 +166.723 5.45125 +166.733 5.7501 +170.021 6.04877 +170.031 6.34745 +170.04 6.64613 +170.05 6.94481 +170.06 7.24349 +170.069 7.54216 +170.079 7.84101 +170.089 8.13969 +173.377 8.43837 +173.387 8.73704 +173.397 9.03572 +173.406 9.3344 +173.416 9.63308 +173.426 9.93192 +173.435 10.2306 +173.445 10.5293 +173.455 10.828 +176.743 11.1266 +176.753 11.4253 +176.763 11.7242 +176.772 12.0228 +176.782 12.3215 +176.792 12.6202 +176.802 12.9189 +176.811 13.2175 +180.1 13.5162 +180.109 13.8151 +180.119 14.1138 +180.129 14.4124 +180.138 14.7111 +180.148 15.0098 +180.158 15.3085 +180.168 15.6071 +180.177 15.906 +180.187 16.2045 +180.197 16.5033 +183.485 16.8019 +183.495 17.1007 +183.504 17.3994 +183.514 17.6982 +183.524 17.9969 +183.534 18.2956 +183.543 18.5943 +183.553 18.8929 +183.563 19.1916 +183.572 19.4903 +183.582 19.7891 +183.592 20.0878 +183.601 20.3865 +183.611 20.6852 +183.621 20.9838 +186.909 21.2825 +186.919 21.5812 +186.929 21.88 +186.938 22.1786 +186.948 22.4774 +186.958 22.7759 +186.967 23.0748 +186.977 23.3733 +186.987 23.6721 +186.997 23.971 +187.006 24.2696 +187.016 24.5683 +190.304 24.867 +190.314 25.1657 +190.324 25.4644 +190.333 25.7632 +190.343 26.0619 +190.353 26.3606 +190.363 26.6592 +190.372 26.9579 +190.382 27.2566 +190.392 27.5553 +190.401 27.8541 +190.411 28.1528 +193.699 28.4515 +193.709 28.7501 +193.719 29.0488 +193.728 29.3475 +193.738 29.6462 +193.748 29.945 +193.758 30.2437 +193.767 30.5424 +193.777 30.8411 +193.787 31.1397 +193.796 31.4384 +197.085 31.7371 +197.095 32.0359 +197.104 32.3346 +197.114 32.6333 +197.124 32.932 +197.133 33.2306 +197.143 33.5293 +197.153 33.8282 +197.162 34.1269 +197.172 34.4255 +197.182 34.724 +197.192 35.0229 +200 35 +200 35 +201 35 +201 35 +202 35 +202 35 +203 35 +203 35 +204 35 +204 35 +205 35 +205 35 +206 35 +206 35 +207 35 +207 35 +208 35 +208 35 +209 35 +209 35 +210 35 +210 35 +211 35 +211 35 +212 35 +212 35 +213 35 +213 35 +214 35 +214 35 +215 35 +215 35 +216 35 +216 35 +217 35 +217 35 +218 35 +218 35 +219 35 +219 35 +220 35 +220 35 +221 35 +221 35 +222 35 +222 35 +223 35 +223 35 +224 35 +224 35 +225 35 +225 35 +226 35 +226 35 +227 35 +227 35 +228 35 +228 35 +229 35 +229 35 +230 35 +230 35 +231 35 +231 35 +232 35 +232 35 +233 35 +233 35 +234 35 +234 35 +235 35 +235 35 +236 35 +236 35 +237 35 +237 35 +238 35 +238 35 +239 35 +239 35 +240 35 +240 35 +241 35 +241 35 +242 35 +242 35 +243 35 +243 35 +244 35 +244 35 +245 35 +245 35 +246 35 +246 35 +247 35 +247 35 +248 35 +248 35 +249 35 +249 35 +250 35 +250 35 +251 35 +251 35 +252 35 +252 35 +253 35 +253 35 +254 35 +254 35 +255 35 +255 35 +256 35 +256 35 +257 35 +257 35 +258 35 +258 35 +259 35 +259 35 +260 35 +260 35 +261 35 +261 35 +262 35 +262 35 +263 35 +263 35 +264 35 +264 35 +265 35 +265 35 +266 35 +266 35 +267 35 +267 35 +268 35 +268 35 +269 35 +269 35 +270 35 +270 35 +271 35 +271 35 +272 35 +272 35 +273 35 +273 35 +274 35 +274 35 +275 35 +275 35 +276 35 +276 35 +277 35 +277 35 +278 35 +278 35 +279 35 +279 35 +280 35 +280 35 +281 35 +281 35 +282 35 +282 35 +283 35 +283 35 +284 35 +284 35 +285 35 +285 35 +286 35 +286 35 +287 35 +287 35 +288 35 +288 35 +289 35 +289 35 +290 35 +290 35 +291 35 +291 35 +292 35 +292 35 +293 35 +293 35 +294 35 +294 35 +295 35 +295 35 +296 35 +296 35 +297 35 +297 35 +298 35 +298 35 +299 35 +299 35 +300 35 +301 34.9 +302 34.8 +303 34.7 +304 34.6 +305 34.5 +306 34.4 +307 34.3 +308 34.2 +309 34.1 +310 34 +311 33.9 +312 33.8 +313 33.7 +314 33.6 +315 33.5 +316 33.4 +317 33.3 +318 33.2 +319 33.1 +320 33 +321 32.9 +322 32.8 +323 32.7 +324 32.6 +325 32.5 +326 32.4 +327 32.3 +328 32.2 +329 32.1 +330 32 +331 31.9 +332 31.8 +333 31.7 +334 31.6 +335 31.5 +336 31.4 +337 31.3 +338 31.2 +339 31.1 +340 31 +341 30.9 +342 30.8 +343 30.7 +344 30.6 +345 30.5 +346 30.4 +347 30.3 +348 30.2 +349 30.1 +350 30 +352.452 29.8669 +355.74 30.0444 +359.019 30.0444 +362.297 30.0444 +365.567 29.8669 +368.845 29.8669 +372.124 29.8669 +375.393 29.6893 +378.671 29.6893 +381.94 29.5118 +385.209 29.3343 +388.478 29.1568 +391.757 29.1568 +395.026 28.9792 +395.017 28.8018 +398.285 28.6242 +401.555 28.4466 +404.824 28.2692 +404.814 28.0917 +408.083 27.9142 +411.352 27.7367 +411.342 27.5592 +411.332 27.3817 +411.323 27.2041 +414.592 27.0266 +417.861 26.8491 +417.851 26.6716 +421.12 26.4941 +421.11 26.3166 +424.379 26.1391 +424.37 25.9615 +427.639 25.784 +427.629 25.6065 +430.898 25.429 +430.888 25.2515 +434.157 25.074 +434.147 24.8963 +434.138 24.7189 +437.407 24.5413 +437.397 24.3639 +440.666 24.1864 +440.656 24.0089 +443.925 23.8314 +443.916 23.6538 +443.906 23.4763 +447.175 23.2988 +447.165 23.1213 +447.155 22.9438 +450.424 22.7663 +450.415 22.5888 +450.405 22.4112 +453.674 22.2337 +453.664 22.0562 +453.654 21.8787 +456.924 21.7012 +456.914 21.5237 +456.904 21.3462 +460.173 21.1686 +460.163 20.991 +460.154 20.8136 +463.423 20.636 +463.413 20.4586 +466.682 20.2811 +466.672 20.1036 +466.663 19.926 +466.653 19.7485 +469.922 19.571 +469.912 19.3935 +469.902 19.216 +473.171 19.0385 +473.162 18.8609 +473.152 18.6834 +476.421 18.5059 +476.411 18.3284 +476.402 18.1509 +479.671 17.9734 +479.661 17.7959 +479.651 17.6183 +482.92 17.4408 +482.911 17.2633 +482.901 17.0858 +486.17 16.9083 +486.16 16.7308 +486.151 16.5533 +489.419 16.3757 +489.41 16.1982 +489.4 16.0207 +492.669 15.8432 +492.659 15.6657 +492.65 15.4882 +492.64 15.3107 +495.909 15.1331 +495.899 14.9556 +495.889 14.7781 +499.159 14.6006 +499.149 14.4231 +499.139 14.2456 +502.408 14.068 +502.398 13.8904 +502.389 13.713 +505.658 13.5354 +505.648 13.358 +505.638 13.1804 +508.907 13.003 +508.898 12.8254 +508.888 12.6479 +512.157 12.4704 +512.147 12.2929 +512.138 12.1154 +515.406 11.9379 +515.397 11.7604 +518.666 11.5828 +518.656 11.4053 +518.646 11.2278 +518.637 11.0503 +521.906 10.8728 +521.896 10.6953 +521.886 10.5178 +525.155 10.3402 +525.145 10.1627 +525.136 9.9851 +528.405 9.8077 +528.395 9.6301 +528.385 9.4527 +528.376 9.2751 +531.645 9.0976 +531.635 8.9201 +531.625 8.7426 +534.894 8.5651 +534.885 8.3876 +534.875 8.2101 +534.865 8.0325 +538.134 7.855 +538.125 7.6775 +538.115 7.5 +541.384 7.3225 +541.374 7.145 +541.364 6.9675 +541.355 6.7899 +544.624 6.6124 +544.614 6.4349 +544.604 6.2574 +547.873 6.0799 +547.864 5.9024 +547.854 5.7248 +551.123 5.5473 +551.113 5.3698 +551.103 5.1923 +551.094 5.0148 +554.363 4.8373 +554.353 4.6598 +557.622 4.4822 +557.612 4.3047 +557.603 4.1272 +560.872 3.9497 +560.862 3.7722 +560.852 3.5947 +564.121 3.4172 +564.111 3.2396 +567.38 3.0621 +567.371 2.8846 +570.64 2.7071 +573.909 2.5295 +573.899 2.3521 +577.168 2.1746 +577.158 1.9969 +580.427 1.8195 +583.696 1.642 +586.965 1.4645 +586.956 1.287 +590.225 1.1095 +593.494 0.932 +596.772 0.932 +600.051 0.932 +603.32 0.7544 +606.599 0.7544 +606.589 0.5769 +609.868 0.5769 +613.137 0.3994 +616.415 0.3994 +619.684 0.2219 diff --git a/qsim.cc b/qsim.cc index bfb83d0..763e99e 100644 --- a/qsim.cc +++ b/qsim.cc @@ -89,7 +89,14 @@ int main(int argc, char** argv){ G4VModularPhysicsList* physlist = factory.GetReferencePhysList("FTFP_BERT"); physlist->SetVerboseLevel(verbose); runManager->SetUserInitialization(physlist); - physlist->RegisterPhysics( new qsimOpticalPhysics() ); + + G4VPhysicsConstructor * opt_phys = new qsimOpticalPhysics(true); + //qsimOpticalPhysics* opt_phys = new qsimOpticalPhysics(true); + physlist->RegisterPhysics((qsimOpticalPhysics *)opt_phys ); + rmmess->SetOptPhys(((qsimOpticalPhysics *)opt_phys)); + + //physlist->RegisterPhysics( new qsimOpticalPhysics() ); + //------------------------------- // UserAction classes diff --git a/src/qsimDetectorConstruction.cc b/src/qsimDetectorConstruction.cc index cad42f1..205c7fd 100644 --- a/src/qsimDetectorConstruction.cc +++ b/src/qsimDetectorConstruction.cc @@ -1,8 +1,12 @@ +//From brads_qsim branch #include "qsimDetectorConstruction.hh" #include "G4SystemOfUnits.hh" +#include "G4PhysicalConstants.hh" + #include "qsimDetector.hh" #include "qsimScintDetector.hh" #include "G4SDManager.hh" + #include "G4Material.hh" #include "G4MaterialTable.hh" #include "G4Element.hh" @@ -10,7 +14,13 @@ #include "G4LogicalBorderSurface.hh" #include "G4LogicalSkinSurface.hh" #include "G4Box.hh" +#include "G4ExtrudedSolid.hh" +#include "G4SubtractionSolid.hh" +#include "G4IntersectionSolid.hh" +#include "G4TwoVector.hh" +#include "G4Polycone.hh" #include "G4Trap.hh" +#include "G4GenericTrap.hh" #include "G4Cons.hh" #include "G4Tubs.hh" #include "G4Trd.hh" @@ -19,51 +29,61 @@ #include "G4ThreeVector.hh" #include "G4Transform3D.hh" #include "G4PVPlacement.hh" -#include "G4OpticalSurface.hh" -#include "G4SubtractionSolid.hh" + +#include +#include +#include +#include + #include "G4VisAttributes.hh" +#include "G4OpticalSurface.hh" +#include "G4NistManager.hh" //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... - void qsimDetectorConstruction::DetModeSet(G4int detMode = 1) { - - fDetMode = detMode; - // 0 is PREX-I design - // 1 is PREX-II prototype (so-called "design 3") - + + fDetMode = detMode; + // 0 is PREX-I design + // 1 is PREX-II prototype (so-called "design 3") + +} + +void qsimDetectorConstruction::LGReflectivitySet(G4double ref = 0.9) { + + fLGReflectivity = ref; } void qsimDetectorConstruction::StandModeSet(G4int standMode = 0) { - + fStandMode = standMode; // 1 cosmic setup (detector, lead, scintillators) // 0 beam setup (detector only) - + } +qsimDetectorConstruction::qsimDetectorConstruction() +{ + + det_x = det_y = det_z = 275*cm; + quartz_x = 1.75*cm; + quartz_y = 7.*cm; + //Change quartz thickness here. + quartz_z = 0.5*cm; + + quartz_zPos = -.0*cm; -qsimDetectorConstruction::qsimDetectorConstruction() { - DetModeSet(); StandModeSet(); +//Position scan + fQuartzPolish = 0.97; + fLGReflectivity = 0.9; + //fDetAngle = 0.*deg; + fDetPosX = 0.*cm; + fDetPosY = 0.*cm; - fQuartzPolish = 0.97; - fDetAngle = 0.*deg; - fDetPosX = 0.*cm; - fDetPosY = 0.*cm; - // fNewStand = false; // Default setting is for the setup to reflect to old cosmic stand. True will go to the new design. Messenger has commands to switch between these at command line or at batch mode run time as well. // fAccBeamStand = false; // Only affects stand components: true deletes the lead block. - - quartz_x = 1.75*cm; // CSC measures in SolidWorks 0.689 x 2.952 x 0.197 cm - quartz_y = 7.5*cm; //2.5 - //Half cm - quartz_z = 0.5*cm; - //One cm - // quartz_z = 0.5*cm; - - quartz_zPos = 0.0*cm; //-.9*cm;//-1.1*cm; //-.9*cm; //-.6*cm; - + cone_rmin1 = 2.1*cm; cone_rmax1 = cone_rmin1+.05*cm; cone_rmin2 = 2.5*cm; // normally 2.5*cm; @@ -71,450 +91,1561 @@ qsimDetectorConstruction::qsimDetectorConstruction() { cone_z = quartz_y+.5*cm; //3 cone_sphi = 0.; cone_fphi = 2*3.1415; - + rin = cone_rmin2; // normally 2.5*cm; rout = rin+.05*cm; lngth = 1.6*cm; // PMT dist. = 2*lngth +1cm (10.4 == 4.5, 6.8 == 2.9) + //Edited 12/16 } + //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... qsimDetectorConstruction::~qsimDetectorConstruction(){;} //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -G4VPhysicalVolume* qsimDetectorConstruction::Construct() { - - // Define materials - - G4double a, z, density; - G4int nelements; - - // Air - G4Element* N = new G4Element("Nitrogen", "N", z=7 , a=14.01*g/mole); - G4Element* O = new G4Element("Oxygen" , "O", z=8 , a=16.00*g/mole); - - G4Material* Air = new G4Material("Air", density=1.29*mg/cm3, nelements=2); - Air->AddElement(N, 70.*perCent); - Air->AddElement(O, 30.*perCent); - - // Quartz - G4Element* Si = new G4Element("Silicon", "Si", z=14 , a=28*g/mole); - - G4Material* Quartz = new G4Material("Quartz", density= 2.203*g/cm3, nelements=2); - Quartz->AddElement(Si, 1); - Quartz->AddElement(O, 2); - - // Aluminum for mirror and stand (need separate materials so that mirror can be reflective) - G4Element* Al = new G4Element("Aluminum", "Al", z=13 , a=27*g/mole); - G4Material* Alu_Mat = new G4Material("Alu_Mat", 2.7*g/cm3, nelements=1); - Alu_Mat->AddElement(Al, 1); - G4Material* Mirror = new G4Material("Mirror", density= 2.7*g/cm3, nelements=1); - Mirror->AddElement(Al, 1); - - // Lead - G4Element* Pb = new G4Element("Lead", "Pb", z=82 , a=207.2*g/mole); - G4Material* Pb_Mat = new G4Material("Pb_Mat", 11.34*g/cm3, nelements=1); - Pb_Mat->AddElement(Pb, 1); - - // Let us make cathode from a special metal (reflectivity 0, efficiency of photoelectrons 25%) - G4Material* CATH = new G4Material("CATH", density= 2.7*g/cm3, nelements=1); - CATH->AddElement(Al, 1); - - - // Define optical property tables - - const G4int nEntries = 205; - - // Array of photon energies - G4double PhotonEnergy[nEntries] = - { 2.4,2.42,2.44,2.46,2.48,2.5,2.52,2.54,2.56,2.58, - 2.6,2.62,2.64,2.66,2.68,2.7,2.72,2.74,2.76,2.78, - 2.8,2.82,2.84,2.86,2.88,2.9,2.92,2.94,2.96,2.98, - 3,3.02,3.04,3.06,3.08,3.1,3.12,3.14,3.16,3.18, - 3.2,3.22,3.24,3.26,3.28,3.3,3.32,3.34,3.36,3.38, - 3.4,3.42,3.44,3.46,3.48,3.5,3.52,3.54,3.56,3.58, - 3.6,3.62,3.64,3.66,3.68,3.7,3.72,3.74,3.76,3.78, - 3.8,3.82,3.84,3.86,3.88,3.9,3.92,3.94,3.96,3.98, - 4,4.02,4.04,4.06,4.08,4.1,4.12 ,4.14,4.16,4.18, //Glass cuts off above 4.135eV, 87 entries - 4.2,4.22,4.24,4.26,4.28,4.3,4.32,4.34,4.36,4.38, - 4.4,4.42,4.44,4.46,4.48,4.5,4.52,4.54,4.56,4.58, - 4.6,4.62,4.64,4.66,4.68,4.7,4.72,4.74,4.76,4.78, - 4.8,4.82,4.84,4.86,4.88,4.9,4.92,4.94,4.96,4.98, // Cut off -> 4.96eV ~ 250nm - 5,5.02,5.04 , 5.06,5.08,5.1,5.12,5.14,5.16,5.18, // 5.04eV = 246 nm is the 30% cutoff, 133 entries - 5.2,5.22,5.24,5.26,5.28,5.3,5.32,5.34,5.36,5.38, - 5.4,5.42,5.44,5.46,5.48,5.5,5.52,5.54,5.56,5.58, - 5.6,5.62,5.64,5.66,5.68,5.7,5.72,5.74,5.76,5.78, - 5.8,5.82,5.84,5.86,5.88,5.9,5.92,5.94,5.96,5.98, - 6,6.02,6.04,6.06,6.08,6.1,6.12,6.14,6.16,6.18, - 6.21,6.29,6.38,6.48,6.57,6.67,6.78,6.87,6.98, - 7.08,7.20,7.32,7.44,7.56,7.69 - }; - - // Cathode quantum efficiency - // Response obtained from the plot of the quantum efficiency as a function of wavelength and then changed to eV for the Bialkali photocathode (synthetic silica) +G4VPhysicalVolume* qsimDetectorConstruction::Construct() +{ + +//File + myfile.open("UVS_45total.txt"); +// outfile.open("points.txt"); +// ------------- Materials ------------- + G4double a, z, density; + G4int nelements; + +// Gasses +// + G4Element* N = new G4Element("Nitrogen", "N", z=7 , a=14.01*g/mole); + G4Element* O = new G4Element("Oxygen" , "O", z=8 , a=16.00*g/mole); + G4Element* elAr = new G4Element("Argon" , "Ar", z=18, a=39.95*g/mole); + G4Element* C = new G4Element("Carbon" , "C", z=6 , a=12.01*g/mole); + G4Element* H = new G4Element("Hydrogen", "H", z=1 , a=1.00*g/mole); + +// G4Material* Air = new G4Material("Air", density=1.29*mg/cm3, nelements=1); +// Air->AddElement(N, 100.*perCent); +// Air->AddElement(O, 30.*perCent); + G4Material* Gas; + G4NistManager* man = G4NistManager::Instance(); + G4Material* N2 = new G4Material("N2", density=.001251*g/cm3, nelements=1); + N2->AddElement(N, 2); + G4Material* Ar = new G4Material("Ar", density=.001662*g/cm3, nelements=1); + Ar->AddElement(elAr, 1); + G4Material* CO2 = man->FindOrBuildMaterial("G4_CARBON_DIOXIDE"); + G4Material* Air = man->FindOrBuildMaterial("G4_AIR"); + +// Quartz +// + G4Element* Si = new G4Element("Silicon", "Si", z=14 , a=28*g/mole); + + G4Material* Quartz = new G4Material("Quartz", density= 2.203*g/cm3, nelements=2); + Quartz->AddElement(Si, 1); + Quartz->AddElement(O, 2); + +// Mirror +// + G4Element* Al = new G4Element("Aluminum", "Al", z=13 , a=27*g/mole); + + G4Element* Pb = new G4Element("Lead", "Pb", z=82 , a=207.2*g/mole); + + G4Material* Alu_Mat = new G4Material("Alu_Mat", 2.7*g/cm3, nelements=1); + Alu_Mat->AddElement(Al, 1); + + G4Material* Pb_Mat = new G4Material("Pb_Mat", 11.34*g/cm3, nelements=1); + Pb_Mat->AddElement(Pb, 1); + + //G4Material* Pb_Mat=Air; // To remove lead bricks, uncomment. + + G4Material* Mirror = new G4Material("Mirror", density= 2.7*g/cm3, nelements=1); + Mirror->AddElement(Al, 1); + + //---Mylar + G4Material* Mylar = new G4Material("Mylar",density = 1.38*g/cm3, nelements = 3); + Mylar->AddElement(C,10); + Mylar->AddElement(H,8); + Mylar->AddElement(O,4); + + //---Aluminized-Mylar + G4Material* AlMylar = new G4Material("AlMylar", density = 0.999999*g/cm3, nelements=4); + AlMylar->AddElement(C,10); + AlMylar->AddElement(H,8); + AlMylar->AddElement(O,4); + AlMylar->AddElement(Al,3); + + //---Paper + G4Material* Paper = new G4Material("Paper", density = 1.5*g/cm3, nelements = 3); + Paper->AddElement(C,6); + Paper->AddElement(H,10); + Paper->AddElement(O,5); + +// Let us make cathode from a special metal (reflectivity 0, efficiency of photoelectrons 25%) + G4Material* CATH = new G4Material("CATH", density= 2.7*g/cm3, nelements=1); + CATH->AddElement(Al, 1); + + +// +// ------------ Generate & Add Material Properties Table ------------ +// + + +const G4int nEntries = 255; + + +//This was here before, but don't know what file qe.txt is. +//Also running macros/vis.mac while setting /gun/particle e- +//would give an error "Photon Energy <= 0" + +G4double PhotonEnergy[nEntries]; +//std::cout<<"trying to open file"; +myfile2.open("qe.txt"); + +//G4double PhotonEnergy[nEntries]; +G4double Efficiency4[nEntries]; +//G4double reflectivity[nEntries]; + +//std::cout<<"before"; + if(myfile2.is_open()) +{ +//std::cout<<"start\n"; +char delim2 = ' '; +std::string inen; +int j = 0; +while(!myfile2.peek() != EOF && j 4.96eV ~ 250nm + 5,5.02,5.04,5.06,5.08,5.1,5.12,5.14,5.16,5.18, // 5.04eV = 246 nm is the 30% cutoff, 133 entries + 5.2,5.22,5.24,5.26,5.28,5.3,5.32,5.34,5.36,5.38, + 5.4,5.42,5.44,5.46,5.48,5.5,5.52,5.54,5.56,5.58, + 5.6,5.62,5.64,5.66,5.68,5.7,5.72,5.74,5.76,5.78, + 5.8,5.82,5.84,5.86,5.88,5.9,5.92,5.94,5.96,5.98, + 6,6.02,6.04,6.06,6.08,6.1,6.12,6.14,6.16,6.18, //200 nm + 6.22, 6.26, 6.3, 6.34, 6.38, 6.42, 6.46, 6.5, + 6.54, 6.58, 6.62, 6.66, 6.7, 6.74, 6.78, 6.82, + 6.86, 6.9, 6.94, 6.98, 7.02, 7.06, 7.1, 7.14, + 7.18, 7.22, 7.26, 7.3, 7.34, 7.38, 7.42, 7.46, + 7.5, 7.54, 7.58, 7.62, 7.66, 7.7, 7.74, 7.78}; */ + + G4double RefractiveIndex1[nEntries]; //quartz ref index + G4double Absorption1[nEntries]; //quartz absorbtion + G4double RefractiveIndex2[nEntries]; //air ref index + G4double RefractiveIndex3[nEntries]; //mirror ref index + G4double Reflectivity4[nEntries]; //cathode reflectivity + G4double Reflectivity3[nEntries]; //mirror reflectivity + G4double ReflectivityPMT[nEntries]; //pmt window reflectivity + G4double RefractiveIndexAR[nEntries]; //argon ref index + G4double RefractiveIndexN2[nEntries]; //nitrogen ref index + G4double RefractiveIndexCO2[nEntries]; //carbon dioxide ref index + + //reflectivity for light guide + G4double Reflect_LG[nEntries]; - G4double EfficiencyArrayPercent[nEntries] = - { 11.0,12.0,12.5,13.1,13.5,14.5,15.2,16.0,16.5,17.0, // percentages here - 17.5,18.0,18.5,19.0,19.2,19.7,20.1,20.7,20.9,21.1, - 21.6,22.0,22.5,22.7,23.0,23.5,23.7,24.0,24.0,24.2, - 24.2,24.5,25.0,25.0,25.3,25.5,25.5,25.5,25.5,25.5, - 25.5,25.5,25.5,25.7,26.1,26.1,26.1,26.1,25.6,25.6, - 25.6,25.6,25.6,25.6,25.6,25.6,25.6,25.6,25.6,25.6, - 25.6,25.6,25.6,26.1,26.1,26.1,26.1,26.1,26.1,26.1, - 25.6,25.0,25.0,25.0,25.0,24.5,24.5,24.5,24.5,24.3, - 24.0,24.0,24.0,24.0,24.0,24.0,23.5,23.5,23.5,23.5, - 23.5,23.5,23.5,23.3,23.1,22.8,22.6,22.6,22.6,22.6, // 4.38 eV - 22.6,22.6,22.3,22.1,22.1,22.1,22.0,21.8,21.7,21.3, //100 entries 4.58 eV - 21.2,21.0,20.8,20.8,20.8,20.8,20.8,20.8,20.8,20.8, // 4.78 eV - 20.4,20.4,20.4,20.4,20.4,20.4,20.4,20.4,20.2,20.0, // 4.98 eV - 20.0,20.0,20.0,20.0,20.0,20.0,20.0,19.5,19.5,19.5, // 5.18 eV - 19.5,19.5,19.5,19.5,19.1,19.1,19.1,19.1,19.1,19.1, // 5.38 eV - 19.1,19.1,19.1,19.0,18.8,18.8,18.8,18.8,18.8,18.8, // 5.58 eV - 18.8,18.4,18.4,18.4,18.4,18.4,18.4,18.4,18.4,18.4, // 5.78 eV - 18.4,18.4,18.4,18.4,18.4,18.4,18.4,18.4,18.4,18.4, // 5.98eV - 18.4,18.2,18.0,18.0,18.0,18.0,18.0,18.0,18.0,18.0, - 18.0,17.6,17.6,17.6,17.6,17.2,16.5,16.2,15.9, - 15.2,14.9,14.3,12.1,10.2,9.6 - }; // 6.18 eV + //reflectivity for blackening + G4double Reflect_Blk[nEntries]; - G4double EfficiencyArray[nEntries]; + //reflectivity paper + G4double Reflect_Paper[nEntries]; - G4double RefractiveIndex1[nEntries]; // quartz refractive index - G4double Absorption1[nEntries]; // quartz absorption length - G4double RefractiveIndex2[nEntries]; // air refractive index - G4double Reflectivity1[nEntries]; // mirror reflectivity - G4double Reflectivity2[nEntries]; // cathode reflectivity + /*G4double Efficiency4[nEntries] = //cathode efficiency + { + 8.7279, 9.047, 9.3778, 9.7207, 10.076, 10.4444, 10.9285, 10.8284, 11.3268, 11.74, + 12.1684, 11.2264, 11.6368, 12.0623, 12.5032, 12.6118, 12.5057, 13.0714, + 13.5483, 14.0427, 12.9653, 13.4394, 13.9306, 14.042, 13.9334, 14.5537, + 15.0847, 15.6351, 14.4455, 14.9737, 15.5211, 16.2048, 16.0916, 16.2041, + 16.0947, 16.7945, 16.6863, 16.7936, 16.6895, 16.7929, 16.6926, 16.7921, + 16.6959, 17.404, 17.3096, 17.4032, 17.3129, 17.4022, 17.3162, 18.0365, + 17.9527, 18.0356, 17.9561, 18.6928, 18.6162, 18.6919, 18.6197, 18.6909, + 18.6233, 19.3721, 19.3078, 19.3712, 19.3115, 19.3702, 19.3152, 19.3693, 19.3189, + 19.3684, 19.3226, 20.0742, 20.0328, 20.0732, 20.0722, 20.0405, 20.0713, 20.0443, + 20.8027, 20.7811, 20.8016, 20.785, 20.8007, 20.789, 20.7997, 20.793, 21.5575, + 21.5573, 21.5566, 21.5554, 21.5655, 21.5545, 21.5695, 21.5535, 22.3389, + 22.3667, 22.3378, 22.371, 22.3367, 22.3753, 22.3357, 22.3346, 22.3838, + 22.3881, 22.3325, 23.1463, 23.2155, 23.1452, 23.2199, 23.1441, 23.2243, + 23.143, 23.1419, 23.2332, 23.2376, 23.1396, 23.1385, 23.2465, 23.1374, + 24.101, 23.9805, 23.9794, 23.9783, 24.1148, 24.1194, 23.976, 23.9748, + 24.1286, 23.9737, 24.1332, 24.1378, 23.9714, 24.1424, 23.9703, 24.147, + 24.1516, 23.968, 23.9668, 24.1608, 23.9657, 24.1654, 24.17, 23.9634, + 23.9623, 24.1793, 23.9611, 24.1839, 23.3353, 22.3017, 21.5156, 21.7265, + 20.9642, 20.0255, 19.3195, 19.5189, 18.6385, 18.834, 18.8376, 17.9807, + 17.3469, 17.5389, 16.9235, 16.1456, 15.5764, 15.0274, 15.2039, 14.153, + 13.0192, 13.0186, 13.1773, 12.5596, 12.117, 12.2688, 11.6899, 11.8383, + 11.2778, 11.423, 10.8803, 11.0222, 10.1272, 9.7702, 9.9003, 9.4259, + 8.7734, 9.2177, 8.8925, 8.4642, 8.1659, 8.2795, 7.6007, 7.0746, 7.1745, + 6.8252, 6.9228, 6.3529, 6.4443, 6.1289, 6.2182, 5.9129, 6, 5.7045, 5.5036, + 5.7894, 5.5853, 5.3096, 5.1227, 4.9424, 5.3893, 5.1992, 5.0159, 4.7682, + 4.6002, 4.4384, 4.8399, 4.6692, 4.5045, 4.2819, 4.1312, 4.3464, 4.1931, + 3.8453, 3.7099, 4.046, 3.9033, 3.7656, 3.5791, 3.6335, 3.453, 3.3314, + 3.2141, 3.101, 3.3822, 3.263, 3.148, 2.9917, 2.8864, 3.0375, 2.9304, 2.7846, + 2.6865, 2.592, 2.8276, 2.7278, 2.6316, 2.5006, 2.4126, 2.3277, 2.4497, + 2.3633, 2.2456, 2.0903, 2.2804, 2.1224, 1.9456, 1.8771, 1.8111, 1.7473 + }; + */ - - for (int i = 0; i < nEntries; i++) { - - PhotonEnergy[i] = PhotonEnergy[i]*eV; - EfficiencyArray[i] = 0.01*EfficiencyArrayPercent[i]; - RefractiveIndex1[i]= 1.438 + (.01197*PhotonEnergy[i]/eV) - (.001955*PhotonEnergy[i]*PhotonEnergy[i]/eV/eV) + (.0004793*PhotonEnergy[i]*PhotonEnergy[i]*PhotonEnergy[i]/eV/eV/eV); - - // *** need to update this - Absorption1[i] = (exp(4.325)*exp(1.191*PhotonEnergy[i]/eV)*exp(-.213*PhotonEnergy[i]*PhotonEnergy[i]/eV/eV)*exp(-.04086*PhotonEnergy[i]*PhotonEnergy[i]*PhotonEnergy[i]/eV/eV/eV))*m; - if (Absorption1[i] > 25*m) { - Absorption1[i] = 25*m; - } - - // *** need to update this - if (PhotonEnergy[i] < 4.135*eV) { - } else if (PhotonEnergy[i] >= 4.135*eV && PhotonEnergy[i] < 6.203*eV) { - Reflectivity1[i] = .6; // .7 - } else { - Reflectivity1[i] = .5; // .6 +//This loop was commented out in the version from github on the brad_qsim branch +//Looks like this takes the last entry and swaps it with the first +//Swaps element (nEntries-1)-i with element i +//I don't understand why this is here so, I'll leave it commented out I guess... +/* G4cout<<"set data"; + G4double temp; + for(int i = 0; i < nEntries/2; i++){ + temp = PhotonEnergy[nEntries-i-1]; + PhotonEnergy[nEntries-1-i] = PhotonEnergy[i]; + PhotonEnergy[i] = temp; + } +*/ + + G4double wav = 0; + G4double wavelength = 0; + G4double *RefIn1; + if(gasType == 1) RefIn1 = RefractiveIndex2; + else if (gasType == 2) RefIn1 = RefractiveIndexN2; + else if (gasType == 3) RefIn1 = RefractiveIndexCO2; + else RefIn1 = RefractiveIndexAR; + std::string str; + G4double inwav = 0; + G4double inref = 0; + char delim = ' '; +// outfile<<"W-len\t\tRIAir\t\tRIN2\t\tRICO2\t\tRIAr\t\tEnergy\t\tQE\n"; + for (int i = 0; i < nEntries; i++) { + RefractiveIndex1[i]= 1.455 -(.005836*PhotonEnergy[i])+(.003374*PhotonEnergy[i]*PhotonEnergy[i]); + PhotonEnergy[i] = PhotonEnergy[i]*eV; + + +//Aluminum +// Reflectivity3[i] = 0; //.6; + +//Aluminum Real +/* if (PhotonEnergy[i] < 4.135*eV) Reflectivity3[i] = .75; // regularly .75, .7 below .56/.53/.46 tunes to 50 PEs + else if (PhotonEnergy[i] >= 4.135*eV && PhotonEnergy[i] < 6.203*eV) Reflectivity3[i] = .7; + else Reflectivity3[i] = .6; */ // .6 + //Reflectivity3[i] = .6; +//ALZAK +// if (PhotonEnergy[i] < 3.26*eV) { +// Reflectivity3[i]=.93; } +// else { Reflectivity3[i] = 0;} + +// No Mirror +// Reflectivity3[i] = 0; + +// Absorption1[i] = 50.*cm; //Uniform + + Absorption1[i] = (exp(4.325)*exp(1.191*PhotonEnergy[i]/eV)*exp(-.213*PhotonEnergy[i]*PhotonEnergy[i]/(eV*eV))*exp(-.04086*PhotonEnergy[i]*PhotonEnergy[i]*PhotonEnergy[i]/(eV*eV*eV)))*m; + if (Absorption1[i] > 25*m) {Absorption1[i] = 25*m;} + + wavelength = 1.2398*eV/PhotonEnergy[i]; + wav = pow( (wavelength), -2); + RefractiveIndex2[i]=1+(.05792105/(238.0185-wav))+(.00167917/(57.362-wav)); + RefractiveIndex3[i]=0; + Reflectivity4[i]=0; + RefractiveIndexAR[i]=1+(2.50141*pow(10.,-3)/(91.012-wav))+(5.00283*pow(10.,-4)/(87.892-wav))+(5.22343*pow(10.,-2)/(214.02-wav)); + RefractiveIndexN2[i]=1+(6.8552*pow(10.,-5))+(3.243157*pow(10.,-2)/(144-wav)); + RefractiveIndexCO2[i]=1+(6.991*pow(10.,-2)/(166.175-wav))+(1.4472*pow(10.,-3)/(79.609-wav))+(6.42941*pow(10.,-5)/(56.3064-wav))+(5.21306*pow(10.,-5)/(46.0196-wav))+(1.46847*pow(10.,-6)/(.0584738-wav)); +/* if( PhotonEnergy[i] > 2.7*eV ){ + Efficiency4[i]=.25; + } else if ( 2.7*eV >= PhotonEnergy[i] ) { + Efficiency4[i]=.2; + } +*/ Efficiency4[i] = Efficiency4[i]/100.0; + ReflectivityPMT[i] = pow((*(RefIn1 + i)-RefractiveIndex1[i])/(*(RefIn1 + i)+RefractiveIndex1[i]),2); + + //LG + Reflect_LG[i] = fLGReflectivity;//0.9; // Default to 0.9... maybe use Reflectivity3.... + + //Blackening + Reflect_Blk[i]=0; + + //Paper + Reflect_Paper[i] = 0; + //outfile<SetOptStat(0); + //gStyle->SetOptFit(0); + Reflectivity3[i] = inref; + } + else { + Reflectivity3[i] = 0.9; } - - RefractiveIndex2[i] = 1.; - Reflectivity2[i] = 0.; - } - - - // Quartz material property table - - G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable(); - myMPT1->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex1,nEntries); - myMPT1->AddProperty("ABSLENGTH", PhotonEnergy, Absorption1, nEntries); - Quartz->SetMaterialPropertiesTable(myMPT1); - - // Air material properties table - // *** need to add absorption length? - G4MaterialPropertiesTable* myMPT2 = new G4MaterialPropertiesTable(); - myMPT2->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex2, nEntries); - Air->SetMaterialPropertiesTable(myMPT2); - - - G4SDManager* SDman = G4SDManager::GetSDMpointer(); - - - + //PhotonEnergy[i] = PhotonEnergy[i]*eV; + + } + myfile.close(); + //outfile.close(); +//G4cout<<"Closed file"; + +//QUARTZ + + G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable(); + myMPT1->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex1,nEntries); + myMPT1->AddProperty("ABSLENGTH", PhotonEnergy, Absorption1, nEntries); + + Quartz->SetMaterialPropertiesTable(myMPT1); + +// +// Air +// + +/* G4double RefractiveIndex2[nEntries] = + { 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, + 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, + 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, + 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, + 1.00, 1.00, 1.00, 1.00 }; */ + + G4MaterialPropertiesTable* myMPT2 = new G4MaterialPropertiesTable(); + myMPT2->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex2, nEntries); + myMPT2->AddConstProperty("SCINTILLATIONYIELD", 25./MeV); + + + const G4int NUMENTRIES = 2; + G4double Scnt_PP[NUMENTRIES] = { 2.*eV, 6*eV }; + + G4double Scnt_FAST[NUMENTRIES] = { 0.5, 0.5 }; + G4double Scnt_SLOW[NUMENTRIES] = { 0.5, 0.5 }; + + myMPT2->AddProperty("FASTCOMPONENT", Scnt_PP, Scnt_FAST, NUMENTRIES); + myMPT2->AddProperty("SLOWCOMPONENT", Scnt_PP, Scnt_SLOW, NUMENTRIES); + + myMPT2->AddConstProperty("RESOLUTIONSCALE", 2.0); + myMPT2->AddConstProperty("FASTTIMECONSTANT", 1.*ns); + myMPT2->AddConstProperty("SLOWTIMECONSTANT", 10.*ns); + myMPT2->AddConstProperty("YIELDRATIO", 1.0); + + Air->SetMaterialPropertiesTable(myMPT2); + + +// +// Ar +// + + G4MaterialPropertiesTable* myMPTAR = new G4MaterialPropertiesTable(); + myMPTAR->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexAR, nEntries); + myMPTAR->AddConstProperty("SCINTILLATIONYIELD", 510./MeV); + + + + myMPTAR->AddProperty("FASTCOMPONENT", Scnt_PP, Scnt_FAST, NUMENTRIES); + myMPTAR->AddProperty("SLOWCOMPONENT", Scnt_PP, Scnt_SLOW, NUMENTRIES); + + myMPTAR->AddConstProperty("RESOLUTIONSCALE", 2.0); + myMPTAR->AddConstProperty("FASTTIMECONSTANT", 1.*ns); + myMPTAR->AddConstProperty("SLOWTIMECONSTANT", 10.*ns); + myMPTAR->AddConstProperty("YIELDRATIO", 1.0); + + Ar->SetMaterialPropertiesTable(myMPTAR); + + +// +//CO2 +// + + const G4int nEC = 41; + G4double PhotonEC[nEC] = { + 2.04358*eV,2.0664*eV,2.09046*eV,2.14023*eV, + 2.16601*eV,2.20587*eV,2.23327*eV,2.26137*eV, + 2.31972*eV,2.35005*eV,2.38116*eV,2.41313*eV, + 2.44598*eV,2.47968*eV,2.53081*eV,2.58354*eV, + 2.6194*eV,2.69589*eV,2.73515*eV,2.79685*eV, + 2.86139*eV,2.95271*eV,3.04884*eV,3.12665*eV, + 3.2393*eV,3.39218*eV,3.52508*eV,3.66893*eV, + 3.82396*eV,3.99949*eV,4.13281*eV,4.27679*eV, + 4.48244*eV,4.65057*eV,4.89476*eV,5.02774*eV, + 5.16816*eV,5.31437*eV,5.63821*eV,5.90401*eV, + 6.19921*eV}; + G4double CO2_1atm_AbsLen[nEC] = { + 70316.5*m, 66796.2*m, 63314.0*m, 56785.7*m, + 53726.5*m, 49381.2*m, 46640.7*m, 44020.0*m, + 39127.2*m, 36845.7*m, 34671.4*m, 32597.4*m, + 30621.3*m, 28743.4*m, 26154.3*m, 23775.1*m, + 22306.7*m, 19526.3*m, 18263.4*m, 16473.0*m, + 14823.5*m, 12818.8*m, 11053.4*m, 9837.32*m, + 8351.83*m, 6747.67*m, 5648.87*m, 4694.87*m, + 3876.99*m, 3150.27*m, 2706.97*m, 2310.46*m, + 1859.36*m, 1568.2*m, 1237.69*m, 1093.38*m, + 962.586*m, 846.065*m, 643.562*m, 520.072*m, + 133.014*m + }; + + G4MaterialPropertiesTable* myMPTC = new G4MaterialPropertiesTable(); + myMPTC->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexCO2, nEntries); + myMPTC->AddConstProperty("SCINTILLATIONYIELD", 5./MeV); + myMPTC->AddProperty("ABSLENGTH", PhotonEC, CO2_1atm_AbsLen, nEC); + + + myMPTC->AddProperty("FASTCOMPONENT", Scnt_PP, Scnt_FAST, NUMENTRIES); + myMPTC->AddProperty("SLOWCOMPONENT", Scnt_PP, Scnt_SLOW, NUMENTRIES); + + + myMPTC->AddConstProperty("RESOLUTIONSCALE", 2.0); + myMPTC->AddConstProperty("FASTTIMECONSTANT", 1.*ns); + myMPTC->AddConstProperty("SLOWTIMECONSTANT", 10.*ns); + myMPTC->AddConstProperty("YIELDRATIO", 1.0); + + CO2->SetMaterialPropertiesTable(myMPTC); + + +// +//N2 +// + + G4MaterialPropertiesTable* myMPTN = new G4MaterialPropertiesTable(); + myMPTN->AddProperty("RINDEX", PhotonEnergy, RefractiveIndexN2, nEntries); + myMPTN->AddConstProperty("SCINTILLATIONYIELD", 140./MeV); + + + + myMPTN->AddProperty("FASTCOMPONENT", Scnt_PP, Scnt_FAST, NUMENTRIES); + myMPTN->AddProperty("SLOWCOMPONENT", Scnt_PP, Scnt_SLOW, NUMENTRIES); + + myMPTN->AddConstProperty("RESOLUTIONSCALE", 2.0); + myMPTN->AddConstProperty("FASTTIMECONSTANT", 1.*ns); + myMPTN->AddConstProperty("SLOWTIMECONSTANT", 10.*ns); + myMPTN->AddConstProperty("YIELDRATIO", 1.0); + + N2->SetMaterialPropertiesTable(myMPTN); + + //Paper material property table + G4MaterialPropertiesTable* table_material_paper = new G4MaterialPropertiesTable(); + table_material_paper->AddProperty("REFLECTIVITY(P)",PhotonEnergy,Reflect_Paper, nEntries); + Paper->SetMaterialPropertiesTable(table_material_paper); + +// +// Mirror (refractive index = 0) +// + + +/* G4double RefractiveIndex3[nEntries] = + { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, + 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, + 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, + 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, + 0.00, 0.00, 0.00, 0.00 }; */ + + + + +// G4MaterialPropertiesTable* myMPT3 = new G4MaterialPropertiesTable(); +// myMPT3->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex3, nEntries); + // myMPT3->AddProperty("ABSLENGTH", PhotonEnergy, Absorption3, nEntries); + // myMPT3->AddProperty("REFLECTIVITY", PhotonEnergy, Reflectivity3, nEntries); + // myMPT3->AddProperty("EFFICIENCY", PhotonEnergy, Efficiency3, nEntries); + +// Mirror->SetMaterialPropertiesTable(myMPT3); + +// +// CATH +// + +/* + G4double Reflectivity4[nEntries] = + { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, + 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, + 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, + 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, + 0.00, 0.00, 0.00, 0.00 }; + + G4double Efficiency4[nEntries] = + {0.25, 0.25, 0.25, 0.25, 0.25, 0.25, + 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, + 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, + 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, + 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, + 0.25, 0.25 }; */ + + + + + //G4MaterialPropertiesTable* myMPT4 = new G4MaterialPropertiesTable(); + //myMPT4->AddProperty("REFLECTIVITY", PhotonEnergy, Reflectivity4,nEntries); + //myMPT4->AddProperty("EFFICIENCY", PhotonEnergy, Efficiency4, nEntries); + + //CATH->SetMaterialPropertiesTable(myMPT4); + + G4SDManager* SDman = G4SDManager::GetSDMpointer(); + //////////////////////////////////////////// // Defining geometry //////////////////////////////////////////// - + // Volumes - + // World and detector box - - double world_x, world_y, world_z; - double det_x, det_y, det_z; - + + G4double world_x, world_y, world_z; + G4double det_x, det_y, det_z; + world_x = world_y = world_z = 275*cm; - det_x = 15*cm; - det_y = 6*cm; - det_z = 6*cm; - + + det_x = 100*cm; + det_y = 100*cm; + det_z = 100*cm; + + //---world G4Box* world_box = new G4Box("World",world_x,world_y,world_z); - - G4LogicalVolume* world_log - = new G4LogicalVolume(world_box,Air,"World",0,0,0); - + + G4LogicalVolume* world_log = new G4LogicalVolume(world_box,Air,"World_log",0,0,0); world_log->SetVisAttributes(G4VisAttributes::Invisible); - - G4VPhysicalVolume* world_phys - = new G4PVPlacement(0,G4ThreeVector(),world_log,"World",0,false,0); - + + G4VPhysicalVolume* world_phys = new G4PVPlacement(0,G4ThreeVector(),world_log,"World",0,false,0); + + //---det box G4Box* det_box = new G4Box("Detector",det_x,det_y,det_z); - - G4LogicalVolume* det_log - = new G4LogicalVolume(det_box,Air,"Detector_log",0,0,0); - + + G4LogicalVolume* det_log = new G4LogicalVolume(det_box,Air,"Detector_log",0,0,0); det_log->SetVisAttributes(G4VisAttributes::Invisible); - - // First, create solids and logical volumes - - // Quartz (defined for all modes) - - G4double q_yLB = quartz_y - (quartz_z); - - G4Trap* quartz_box = new G4Trap("Quartz", 2*quartz_x, 2*quartz_z, 2*quartz_y, 2*q_yLB); + ///-------------------single cut quartz + ///---Moller style quartz - G4LogicalVolume* quartz_log - = new G4LogicalVolume(quartz_box,Quartz,"Quartz",0,0,0); - - qsimScintDetector* quartzSD = new qsimScintDetector("QuartzSD", 10); - + + //---quartz + G4double q_lz=8.4*cm; //moller quartz + G4double q_ly=2.0*cm; // 1.5cm thick nominal + G4double q_lx_long=16.0*cm; + G4double q_lx_short=q_lx_long-q_ly; // we are doing 45 degree single cut + G4Trap* quartz_box = new G4Trap("Quartz", q_lz, q_ly, q_lx_long, q_lx_short); + G4LogicalVolume* quartz_log = new G4LogicalVolume(quartz_box,Quartz,"Quartz_log",0,0,0); + + //make quartz sensitive detector + qsimScintDetector* quartzSD = new qsimScintDetector("quartzSD", 10); // detector ID is 10 SDman->AddNewDetector(quartzSD); quartz_log->SetSensitiveDetector(quartzSD); - - G4RotationMatrix* rotQ = new G4RotationMatrix; - - rotQ->rotateX(M_PI/2.*rad); - if(fDetMode == 0) { - rotQ->rotateZ(M_PI*rad); - } - - G4VPhysicalVolume* quartz_phys - = new G4PVPlacement(rotQ,G4ThreeVector(0,0,quartz_zPos),quartz_log,"Quartz", det_log,false,0); - // Light guide and tube mirror (only for PREX-I design) - - G4Trap *lightguide_big = new G4Trap("lightguide_big",11.21*cm,2.97456*deg,90.0*deg,3.15*cm+0.05*cm,3.25*cm+0.05*cm,3.25*cm+0.05*cm,0.0*deg,0.225*cm+0.05*cm,2.0*cm+0.05*cm,2.0*cm+0.05*cm,0.0*deg); - G4Trap *lightguide_small = new G4Trap("lightguide_small",11.22*cm,2.97456*deg,90.0*deg,3.15*cm,3.25*cm,3.25*cm,0.0*deg,0.225*cm,2.0*cm,2.0*cm,0.0*deg); - - G4VSolid *lightguide_virt = new G4SubtractionSolid("lightguide_virt", lightguide_big, lightguide_small); - - G4LogicalVolume *lightguide_log = new G4LogicalVolume(lightguide_virt, Mirror, "LG_log",0,0,0); - - G4Cons* mirror_tube = new G4Cons("TMirror",cone_rmin2,cone_rmax2,2.5*cm,2.55*cm,lngth,cone_sphi,cone_fphi); - - G4LogicalVolume* tmirror_log = new G4LogicalVolume(mirror_tube,Mirror,"TMirror",0,0,0); - + //put quartz in det box + G4RotationMatrix* rotate_quartz = new G4RotationMatrix; + rotate_quartz->rotateX(-90.0*M_PI/180.0*rad); // cut part is facing the beam dump + G4VPhysicalVolume* quartz_phys = new G4PVPlacement(rotate_quartz,G4ThreeVector(0,0,0),quartz_log,"Quartz", det_log,false,0); //quartz center is on (0,0,0) - - // PMT and cathode - - G4double anini = 0.*deg; - G4double anspan = 360.*deg; - - G4double prin = 0; - G4double prout = 2.6*cm; - G4double plngth = 1.5*mm; - - G4Tubs* pmt = new G4Tubs("PMT",prin,prout,plngth,anini,anspan); - - G4LogicalVolume* pmt_log - = new G4LogicalVolume(pmt,Air,"PMT",0,0,0); - G4String DetSDname = "tracker1"; - - G4double cin = 0; - G4double cout = 2.6*cm; - G4double clngth = 0.1*mm; - - G4Tubs* cath = new G4Tubs("CATH",cin,cout,clngth,anini,anspan); - - G4LogicalVolume* cath_log - = new G4LogicalVolume(cath,CATH,"CATH",0,0,0); - - qsimDetector* cathSD = new qsimDetector("cath", 2); - - SDman->AddNewDetector(cathSD); - cath_log->SetSensitiveDetector(cathSD); - - G4VisAttributes *cathatt = new G4VisAttributes(); - cathatt->SetColour(1.0, 1.0, 0.2); - cath_log->SetVisAttributes(cathatt); - - - // Scintillators - // Coincidence volumes **** NOTE: Upper scint is below the quartz (First coincidence w/ e-) - - - G4Box* upperScint = new G4Box("upperScint",4.5*cm,4.5*cm,0.75*cm); - G4LogicalVolume* uScint_log = new G4LogicalVolume(upperScint,Air,"upperScint",0,0,0); - - // Make sensitive - DetSDname = "tracker2"; - - qsimScintDetector* upScintSD = new qsimScintDetector(DetSDname, 1); - - SDman->AddNewDetector(upScintSD); - uScint_log->SetSensitiveDetector(upScintSD); - G4double upScint_pos; - - upScint_pos = quartz_z-50*cm; //45*cm; // changed to 45 cm from 50 cm as a rough estimate based on CAD measurements of the PMT model 1 + quartz design on the new stand design. - - - - - G4Box* lowScint = new G4Box("lowScint",4.5*cm,4.5*cm,0.75*cm); - G4LogicalVolume* lScint_log = new G4LogicalVolume(lowScint,Air,"lowScint",0,0,0); - - // Make sensitive - DetSDname = "tracker3"; - - qsimScintDetector* loScintSD = new qsimScintDetector(DetSDname, 2); - - SDman->AddNewDetector(loScintSD); - lScint_log->SetSensitiveDetector(loScintSD); + ///-------------------single cut quartz + ///---Super-elastic style quartz +/* - G4double loScint_pos; - loScint_pos = upScint_pos+1.006*m; //new setup is 1.02874*m; // measured to be 1.02874 m between the two scintillators in the CAD drawings. Previously was just 1.0 m + //---quartz + G4double q_lz=17.1*cm; //moller quartz + G4double q_ly=2.0*cm; // 1.5cm thick nominal + G4double q_lx_long=4.0*cm; + G4double q_lx_short=q_lx_long-q_ly; // we are doing 45 degree single cut + G4Trap* quartz_box = new G4Trap("Quartz", q_lz, q_ly, q_lx_long, q_lx_short); + G4LogicalVolume* quartz_log = new G4LogicalVolume(quartz_box,Quartz,"Quartz_log",0,0,0); - - - // LEAD BLOCK - /////////////////////////////////////////////////////////////////////////////////////// - - G4Box* Pb_blox = new G4Box("Pb_blox", 10.16*cm,7.62*cm, 10.16*cm); - // expanded to ensure nothing - // can hit the scint. w/o the lead. - G4LogicalVolume* Pb_log = new G4LogicalVolume(Pb_blox,Pb_Mat,"Lead",0,0,0); - - G4double Pb_pos; - Pb_pos = loScint_pos-15.35*cm; // new setup is = loScint_pos-18.554*cm; //(-1*quartz_z)+(30.0*cm-(quartz_y*sin(scintAngle)))*sin(scintAngle); - // If fAccBeamStand == true then remove the lead bricks, else leave them + //make quartz sensitive detector + qsimScintDetector* quartzSD = new qsimScintDetector("quartzSD", 10); // detector ID is 10 + SDman->AddNewDetector(quartzSD); + quartz_log->SetSensitiveDetector(quartzSD); - // Detector + //put quartz in det box + G4RotationMatrix* rotate_quartz = new G4RotationMatrix; + rotate_quartz->rotateX(-M_PI/2.0*rad); // cut part is facing the beam dump + G4VPhysicalVolume* quartz_phys = new G4PVPlacement(rotate_quartz,G4ThreeVector(0,0,0),quartz_log,"Quartz", det_log,false,0); //quartz center is on (0,0,0) - - - - // Place physical volumes - - G4RotationMatrix* detrot = new G4RotationMatrix; - G4RotationMatrix* rot_pmt = new G4RotationMatrix; - if (fDetMode == 0) { - - detrot->rotateY(45.*deg); - G4RotationMatrix* rotlg = new G4RotationMatrix; - rotlg->rotateY(M_PI/2.*rad); - rotlg->rotateZ(-M_PI/2.*rad); - - rot_pmt->rotateY(M_PI/2.*rad); - G4VPhysicalVolume* tmirror_phys = new G4PVPlacement(rot_pmt,G4ThreeVector(7.25*cm+lngth+2.*cm,0.,.9*cm),tmirror_log,"TMirror",det_log,false,0); - G4VPhysicalVolume* lightguide_phys = new G4PVPlacement(rotlg,G4ThreeVector(0.*cm,0,-0.375*cm+.9*cm),lightguide_log,"lightguide_phys", det_log,false,0); - G4VPhysicalVolume* pmt_phys = new G4PVPlacement(rot_pmt,G4ThreeVector(7.25*cm+plngth+7.*cm,0.,.9*cm),pmt_log,"PMT",det_log,false,0); - G4VPhysicalVolume* cath_phys = new G4PVPlacement(rot_pmt,G4ThreeVector(7.25*cm+2.*plngth+7.*cm,0.,.9*cm),cath_log,"CATH",det_log,false,0); - } +*/ - if (fDetMode == 1) { - - detrot->rotateY(fDetAngle); - rot_pmt -> rotateY(M_PI/4.*rad); - G4VPhysicalVolume* pmt_phys = new G4PVPlacement(rot_pmt,G4ThreeVector(7.5*cm,0.*cm,0.*mm),pmt_log,"PMT",det_log,false,0); - G4VPhysicalVolume* cath_phys = new G4PVPlacement(rot_pmt,G4ThreeVector(7.5*cm+plngth*cos(M_PI/4*rad)*mm,0.*cm,-plngth*cos(M_PI/4.*rad)*mm),cath_log,"CATH",det_log,false,0); - } +/* - G4VPhysicalVolume* det_phys - = new G4PVPlacement(detrot,G4ThreeVector(fDetPosX,fDetPosY,0.*cm),det_log,"detector_phys",world_log,false,0); - - - if (fStandMode == 1) { - G4PVPlacement* uScint_phys = new G4PVPlacement(0,G4ThreeVector(0.0,0.0,upScint_pos-1.*cm),uScint_log,"upperScint",world_log,false,0); - G4PVPlacement* lScint_phys = new G4PVPlacement(0,G4ThreeVector(0.0,0.0,loScint_pos),lScint_log,"lowerScint",world_log,false,0); - G4PVPlacement* Pb_phys = new G4PVPlacement(0,G4ThreeVector(0.0,0.0*cm,Pb_pos),Pb_log,"Pb",world_log,false,0); - } + ///--------------double cut quartz + ///---Moller style quartz + + //---quartz + //define the poligon + G4double thickness=1.0*cm; + G4double x1=8.0*cm; + G4double y1=0.0*cm; + + G4double x2=8.0*cm - thickness/2; + G4double y2=thickness/2; + + G4double x3=-8.0*cm; + G4double y3=y2; + + G4double x4=-8.0*cm; + G4double y4=-1.0*y3; + + G4double x5=x2; + G4double y5=-1.0*y2; + + std::vector bar_poligon(5); + bar_poligon[0]= G4TwoVector(x1,y1); + bar_poligon[1]= G4TwoVector(x2,y2); + bar_poligon[2]= G4TwoVector(x3,y3); + bar_poligon[3]= G4TwoVector(x4,y4); + bar_poligon[4]= G4TwoVector(x5,y5); + + G4ExtrudedSolid *quartz_box=new G4ExtrudedSolid("Quartz",bar_poligon, 4.2*cm, G4TwoVector(0,0),1.0, G4TwoVector(0,0), 1.0 ); + + G4LogicalVolume* quartz_log = new G4LogicalVolume(quartz_box,Quartz,"Quartz_log",0,0,0); + + //make quartz sensitive detector + qsimScintDetector* quartzSD = new qsimScintDetector("quartzSD", 10); // detector ID is 10 + SDman->AddNewDetector(quartzSD); + quartz_log->SetSensitiveDetector(quartzSD); + + //put quartz in det box + G4RotationMatrix* rotate_quartz = new G4RotationMatrix; + rotate_quartz->rotateX(-M_PI/2.0*rad); // cut part is facing the beam dump + G4VPhysicalVolume* quartz_phys = new G4PVPlacement(rotate_quartz,G4ThreeVector(0,0,0),quartz_log,"Quartz", det_log,false,0); //quartz center is on (0,0,0) + +*/ + +/* + ///--------------double cut quartz + ///---super elastic style quartz + + //---quartz + //define the poligon + G4double thickness=2.0*cm; + G4double x1=2.0*cm; + G4double y1=0.0*cm; + + G4double x2=2.0*cm - thickness/2; + G4double y2=thickness/2; + + G4double x3=-2.0*cm; + G4double y3=y2; + + G4double x4=-2.0*cm; + G4double y4=-1.0*y3; + + G4double x5=x2; + G4double y5=-1.0*y2; + + std::vector bar_poligon(5); + bar_poligon[0]= G4TwoVector(x1,y1); + bar_poligon[1]= G4TwoVector(x2,y2); + bar_poligon[2]= G4TwoVector(x3,y3); + bar_poligon[3]= G4TwoVector(x4,y4); + bar_poligon[4]= G4TwoVector(x5,y5); + + G4ExtrudedSolid *quartz_box=new G4ExtrudedSolid("Quartz",bar_poligon, 17.2/2.0*cm, G4TwoVector(0,0),1.0, G4TwoVector(0,0), 1.0 ); + + G4LogicalVolume* quartz_log = new G4LogicalVolume(quartz_box,Quartz,"Quartz_log",0,0,0); + + //make quartz sensitive detector + qsimScintDetector* quartzSD = new qsimScintDetector("quartzSD", 10); // detector ID is 10 + SDman->AddNewDetector(quartzSD); + quartz_log->SetSensitiveDetector(quartzSD); + + //put quartz in det box + G4RotationMatrix* rotate_quartz = new G4RotationMatrix; + rotate_quartz->rotateX(-M_PI/2.0*rad); // cut part is facing the beam dump + G4VPhysicalVolume* quartz_phys = new G4PVPlacement(rotate_quartz,G4ThreeVector(0,0,0),quartz_log,"Quartz", det_log,false,0); //quartz center is on (0,0,0) + +*/ + + G4double in = 2.54*cm; + + G4double thick_ness = 0.47625*cm; + + G4double x_1 = 7*cm; + G4double y_1 = -1*cm; //-0.5*cm; + + G4double x_2 = -1*cm; + G4double y_2 = -1*cm; //-0.5*cm; + + G4double x_3 = -1*cm; + G4double y_3 = 1*cm;//0.5*cm; + + G4double x_4 = 7*cm; + G4double y_4 = 1*cm;//0.5*cm; + + G4double x_5 = 7*cm + 8*cos(11.5*M_PI/180*rad)*cm; //7*cm + 8*cos(11.5*M_PI/180*rad)*cm; + G4double y_5 = 1*cm + 8*sin(11.5*M_PI/180*rad)*cm; //0.5*cm + 8*sin(11.5*M_PI/180*rad)*cm; + + G4double x_6 = x_5 + 3.0*cm;//+3.5*cm + G4double y_6 = y_5; + + G4double x_7 = x_6; + G4double y_7 = y_6 -8.4*cm; + + G4double x_8 = x_1 + 9.19*cos(36.4822801*deg)*cm;//x_1 + 9.75*cos(36.4822801*deg)*cm + G4double y_8 = y_1 - 9.19*sin(36.4822801*deg)*cm;//y_1 - 9.75*sin(36.4822801*deg)*cm + + G4double poss = q_lz/2 + thick_ness/2; + + std::vector bar_polygon(8); + bar_polygon[0] = G4TwoVector(x_1,y_1); + bar_polygon[1] = G4TwoVector(x_2,y_2); + bar_polygon[2] = G4TwoVector(x_3,y_3); + bar_polygon[3] = G4TwoVector(x_4,y_4); + bar_polygon[4] = G4TwoVector(x_5,y_5); + bar_polygon[5] = G4TwoVector(x_6,y_6); + bar_polygon[6] = G4TwoVector(x_7,y_7); + bar_polygon[7] = G4TwoVector(x_8,y_8); + + + + G4ExtrudedSolid* cone_side = new G4ExtrudedSolid("Cone_Side",bar_polygon,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); + G4LogicalVolume* cone_log = new G4LogicalVolume(cone_side,AlMylar,"Cone_log",0,0,0); + + G4RotationMatrix* rotcone_side = new G4RotationMatrix(); + rotcone_side->rotateX(-M_PI/2*rad); + G4VPhysicalVolume* cone_phys1 = new G4PVPlacement(rotcone_side,G4ThreeVector(0,poss,0),cone_log,"Cone side1",det_log,false,0); + + + G4VPhysicalVolume* cone_phys2 = new G4PVPlacement(rotcone_side,G4ThreeVector(0,-poss,0),cone_log,"Cone side2",det_log,false,0); + + + G4VisAttributes *con = new G4VisAttributes(); + con->SetColour(1.0,0.0,1.0); + cone_log->SetVisAttributes(con); + +/* +//paper for quartz +G4double px1 = q_lx_short/2; +G4double py1 = q_lz/2; +G4double pz1 = in/1000; + +G4Box* papert = new G4Box("papert",px1,py1,pz1); +G4LogicalVolume* papert_l = new G4LogicalVolume(papert,Paper,"PaperL",0,0,0); +G4VPhysicalVolume* papert_p = new G4PVPlacement(0,G4ThreeVector(0,0,q_ly/2+pz1),papert_l,"Paper Top",det_log,false,0); + +G4VisAttributes *paper = new G4VisAttributes(); +paper->SetColour(0.0,0.0,0.0); +papert_l->SetVisAttributes(paper); + +G4double px2 = q_lx_long/2; +G4double py2 = q_lz/2; +G4double pz2 = in/1000; + +G4Box* paperb = new G4Box("paperb",px2,py2,pz2); +G4LogicalVolume* paperb_l = new G4LogicalVolume(paperb,Paper,"PaperL1",0,0,0,true); +G4VPhysicalVolume* paperb_p = new G4PVPlacement(0,G4ThreeVector(0,0,-q_ly/2),paperb_l,"Paper Bot",det_log,false,0); + +paperb_l->SetVisAttributes(paper); + +G4double px3 = q_lz/2; +G4double py3 = q_ly/2; +G4double pz3 = pz2; + +G4Box* paperba = new G4Box("paperba",px2,py2,pz2); +G4LogicalVolume* paperba_l = new G4LogicalVolume(paperba,Paper,"PaperL1",0,0,0,true); + +G4VPhysicalVolume* paperba_p = new G4PVPlacement(0,G4ThreeVector(-q_lx_short/2,0,0),paperba_l,"Paper Bot",det_log,false,0); + +paperba_l->SetVisAttributes(paper); + +G4Box* papers = new G4Box("papers",px2,py3,pz3); +G4LogicalVolume* papers_l = new G4LogicalVolume(papers,Paper,"Papers",0,0,0,true); +G4VPhysicalVolume* papers_p1 = new G4PVPlacement(rotcone_side,G4ThreeVector(0,q_lz/2,0),papers_l,"PaperSide1",det_log,false,0); +G4VPhysicalVolume* papers_p2 = new G4PVPlacement(rotcone_side,G4ThreeVector(0,-q_lz/2,0),papers_l,"PaperSide2",det_log,false,0); + +papers_l->SetVisAttributes(paper); +*/ + + +//reflector + G4double xone = x_4; + G4double yone = q_lz/2 + thick_ness; + + G4double xtwo = x_5; + G4double ytwo = q_lz/2 + thick_ness; + + G4double xthr = x_5; + G4double ythr = -q_lz/2 - thick_ness; + + G4double xfo = x_4; + G4double yfo = -q_lz/2 - thick_ness; + + + std::vector barpolygon(4); + barpolygon[0] = G4TwoVector(xone,yone); + barpolygon[1] = G4TwoVector(xtwo,ytwo); + barpolygon[2] = G4TwoVector(xthr,ythr); + barpolygon[3] = G4TwoVector(xfo,yfo); + + G4ExtrudedSolid* ref_box = new G4ExtrudedSolid("REF",barpolygon,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); + + + G4LogicalVolume* ref_log = new G4LogicalVolume(ref_box,AlMylar,"Ref_log",0,0,0); + + G4RotationMatrix* rot_ref = new G4RotationMatrix(); + rot_ref->rotateY(11.5*M_PI/180*rad); + + G4VPhysicalVolume* refl_phy = new G4PVPlacement(rot_ref,G4ThreeVector(0,0,-thick_ness),ref_log,"Reflector",det_log,false,0); + + + G4VisAttributes *reflect = new G4VisAttributes(); + reflect->SetColour(0.0,1.0,1.0); + ref_log->SetVisAttributes(reflect); +/* +//cone bottom +G4double xuno = x_1; +G4double yuno = q_lz/2 + thick_ness; + +G4double xdos = x_8; +G4double ydos = q_lz/2 + thick_ness; + +G4double xtres = x_8; +G4double ytres = -q_lz/2 -thick_ness; + +G4double xcuat = x_1; +G4double ycuat = -q_lz/2 -thick_ness; + +std::vector barrpolygon(4); +barrpolygon[0] = G4TwoVector(xuno,yuno); +barrpolygon[1] = G4TwoVector(xdos,ydos); +barrpolygon[2] = G4TwoVector(xtres,ytres); +barrpolygon[3] = G4TwoVector(xcuat,ycuat); + + +G4ExtrudedSolid* cone_bot = new G4ExtrudedSolid("Bottom",barrpolygon,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); +*/ + + G4double xbox = (9.35/2)*cm; + G4double ybox = (9.19/2)*cm; //(9.75/2)*cm; + G4double zbox = thick_ness/2; + + G4Box* cone_bot = new G4Box("bottom",xbox,ybox,zbox); + G4LogicalVolume* coneb_log = new G4LogicalVolume(cone_bot,AlMylar,"Coneb_log",0,0,0,true); + + + G4RotationMatrix* rot_coneb = new G4RotationMatrix(); + rot_coneb->rotateY(-36.4822801*M_PI/180*rad); + + + G4VPhysicalVolume* coneb_phy = new G4PVPlacement(rot_coneb,G4ThreeVector(x_1+ ybox*cos(36.4822801*deg),0,-q_ly-ybox*sin(36.4822801*deg)),coneb_log,"coneback",det_log,false,0); + + + + coneb_log->SetVisAttributes(con); + + + + +//another cone top part + + std::vector bpol1(4); + bpol1[0] = G4TwoVector(x_5,q_lz/2+ thick_ness); + bpol1[1] = G4TwoVector(x_6,q_lz/2+thick_ness); + bpol1[2] = G4TwoVector(x_6,-(q_lz/2+thick_ness)); + bpol1[3] = G4TwoVector(x_5, -(q_lz/2+thick_ness)); + + + G4ExtrudedSolid* con_par1 = new G4ExtrudedSolid("Par",bpol1,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); + G4LogicalVolume* con_log1 = new G4LogicalVolume(con_par1,AlMylar,"Conepar",0,0,0,true); + G4VPhysicalVolume* con_phy11 = new G4PVPlacement(0,G4ThreeVector(0,0,y_5+thick_ness/2),con_log1, "cone part11",det_log,false,0); + G4VPhysicalVolume* con_phy22 = new G4PVPlacement(0,G4ThreeVector(0,0,y_7-thick_ness/2),con_log1, "cone part22",det_log,false,0); + + con_log1->SetVisAttributes(con); + + + +//Light guide +//Light guide side + G4double lx1 = x_6; + G4double ly1 = y_6; + + G4double lx2 = x_6 + 2.94*cm; + G4double ly2 = y_6; + + G4double lx3 = lx2 + 39.43217011394904*cm; + G4double ly3 = y_6 - 8.022568094008*cm; + + G4double lx4 = lx1 + 2*cm + (39.3*cos(11.5*M_PI/180*rad))*cm; + G4double ly4 = ly1 -9.35*cm - (39.3*sin(11.5*M_PI/180*rad))*cm; + + G4double lx5 = lx1 + 2*cm; + G4double ly5 = ly1 - 9.35*cm; + + G4double lx6 = lx1; + G4double ly6 = ly1 - 9.35*cm; + + std::vector poly(6); + poly[0] = G4TwoVector(lx1,ly1); + poly[1] = G4TwoVector(lx2,ly2); + poly[2] = G4TwoVector(lx3,ly3); + poly[3] = G4TwoVector(lx4,ly4); + poly[4] = G4TwoVector(lx5,ly5); + poly[5] = G4TwoVector(lx6,ly6); + + G4ExtrudedSolid* l_g_side = new G4ExtrudedSolid("LGS",poly,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); + G4LogicalVolume* l_g_slog = new G4LogicalVolume(l_g_side,AlMylar,"LGSL",0,0,0,true); + + G4RotationMatrix* rot_l = new G4RotationMatrix(); + rot_l->rotateX(-M_PI/2*rad); + + G4VPhysicalVolume* l_g_phy1 = new G4PVPlacement(rot_l,G4ThreeVector(-(x_6 - x_5),q_lz/2 + 1.5*thick_ness,thick_ness),l_g_slog,"LGS1",det_log,false,0); + G4VPhysicalVolume* l_g_phy2 = new G4PVPlacement(rot_l,G4ThreeVector(-(x_6 - x_5),-q_lz/2 - 1.5*thick_ness,thick_ness),l_g_slog,"LGS2",det_log,false,0); + + G4VisAttributes* lg = new G4VisAttributes(); + lg->SetColour(0.0,0.0,1.0); + l_g_slog->SetVisAttributes(lg); + + +//Light guide top- part 1 + G4double xxuno = lx1; + G4double yyuno = poss + 1.5*thick_ness; + + G4double xxdos = lx2; + G4double yydos = poss + 1.5*thick_ness; + + G4double xxtres = lx2; + G4double yytres = -poss - 1.5*thick_ness; + + G4double xxcuat = lx1; + G4double yycuat = -poss - 1.5*thick_ness; + + + std::vector topp(4); + topp[0] = G4TwoVector(xxuno,yyuno); + topp[1] = G4TwoVector(xxdos,yydos); + topp[2] = G4TwoVector(xxtres,yytres); + topp[3] = G4TwoVector(xxcuat,yycuat); + + + G4ExtrudedSolid* lg_top = new G4ExtrudedSolid("LGT",topp,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); + G4LogicalVolume* lgt_log = new G4LogicalVolume(lg_top,Mylar,"LGTL",0,0,0,true); + G4VPhysicalVolume* ltb_phy = new G4PVPlacement(0,G4ThreeVector(-(x_6-x_5+thick_ness),0,y_5+thick_ness),lgt_log,"LGBP",det_log,false,0); + + + G4VisAttributes* lgt = new G4VisAttributes(); + lgt->SetColour(0.0,0.0,1.0); + lgt_log->SetVisAttributes(lgt); + + +//Light guide top - part 2 + G4double xxun = lx3 +2*thick_ness; + G4double yyun = poss + 1.5*thick_ness; + + G4double xxdeux = lx2; + G4double yydeux = poss + 1.5*thick_ness; + + G4double xxtrois = lx2; + G4double yytrois = -poss - 1.5*thick_ness; + + G4double xxquatre = lx3 + 2*thick_ness; + G4double yyquatre = -poss - 1.5*thick_ness; + + + std::vector topb(4); + topb[0] = G4TwoVector(xxun,yyun); + topb[1] = G4TwoVector(xxdeux,yydeux); + topb[2] = G4TwoVector(xxtrois,yytrois); + topb[3] = G4TwoVector(xxquatre,yyquatre); + + G4ExtrudedSolid* lg_top1 = new G4ExtrudedSolid("LGT",topb,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); + G4LogicalVolume* lgt_log1 = new G4LogicalVolume(lg_top1,AlMylar,"LGTL",0,0,0,true); + + + G4RotationMatrix *rot_ll = new G4RotationMatrix(); + rot_ll->rotateY(-11.5*M_PI/180*rad); + + G4VPhysicalVolume* ltb_phy1 = new G4PVPlacement(rot_ll,G4ThreeVector(-2*cm -1.75*thick_ness,0,poss+ 5.25*thick_ness),lgt_log1,"LGBP",det_log,false,0); + + + + lgt_log1->SetVisAttributes(lgt); + + + +//Light guide bottom- Part 1 + G4double xx1 = lx5; + G4double yy1 = -poss - 1.5*thick_ness; + + G4double xx2 = lx4 + 1.5*thick_ness; + G4double yy2 = -poss - 1.5*thick_ness; + + G4double xx3 = lx4 + 1.5*thick_ness; + G4double yy3 = poss + 1.5*thick_ness; + + G4double xx4 = lx5; + G4double yy4 = poss + 1.5*thick_ness; + + + std::vector botp(4); + botp[0] = G4TwoVector(xx1,yy1); + botp[1] = G4TwoVector(xx2,yy2); + botp[2] = G4TwoVector(xx3,yy3); + botp[3] = G4TwoVector(xx4,yy4); + + G4ExtrudedSolid* lg_bot = new G4ExtrudedSolid("LGB",botp,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); + G4LogicalVolume* lgb_log = new G4LogicalVolume(lg_bot,AlMylar,"LGBL",0,0,0,true); + + G4VPhysicalVolume* lgb_phy = new G4PVPlacement(rot_ll,G4ThreeVector((-2*cm - thick_ness),0,y_7 + 7*thick_ness),lgb_log,"LGB",det_log,false,0); + + + G4VisAttributes* lgb = new G4VisAttributes(); + lgb->SetColour(0.0,0.0,1.0); + lgb_log->SetVisAttributes(lgb); + + + +//Light guide bottom- Part 2 + G4double xx11 = lx5; + G4double yy11 = -poss - 1.5*thick_ness; + + G4double xx22 = lx6; + G4double yy22 = -poss - 1.5*thick_ness; + + G4double xx33 = lx6; + G4double yy33 = poss + 1.5*thick_ness; + + G4double xx44 = lx5; + G4double yy44 = poss + 1.5*thick_ness; + + + std::vector botp1(4); + botp1[0] = G4TwoVector(xx11,yy11); + botp1[1] = G4TwoVector(xx22,yy22); + botp1[2] = G4TwoVector(xx33,yy33); + botp1[3] = G4TwoVector(xx44,yy44); + + G4ExtrudedSolid* lg_bot2 = new G4ExtrudedSolid("LGB",botp1,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); + G4LogicalVolume* lgb_log2 = new G4LogicalVolume(lg_bot2,AlMylar,"LGBL",0,0,0,true); + + G4VPhysicalVolume* lgb_phy2 = new G4PVPlacement(0,G4ThreeVector(-(x_6-x_5),0,-q_lz+ 3*thick_ness),lgb_log2,"LGBP",det_log,false,0); + + + lgb_log2->SetVisAttributes(lgb); + +/* +G4double pos_s = 4.2*cm + thick_ness; + +//Double cut 1 cm Moller cone sides +G4double x_11 = 7.0*cm; +G4double y_11 = -0.5*cm; + +G4double x_22 = -1.0*cm; +G4double y_22 = -0.5*cm; + +G4double x_33 = -1.0*cm; +G4double y_33 = 0.5*cm; + +G4double x_44 = 7.0*cm; +G4double y_44 = 0.5*cm; + +G4double x_55 = x_44 + 8*cos(13.5*M_PI/180*rad)*cm; + +G4double y_55 = y_44 + 8*sin(13.5*M_PI/180*rad)*cm -0.5*cm; + +G4double x_66 = x_55; +G4double y_66 = y_55 + 0.5*cm; + +G4double x_77 = x_66 + 2.0*cm; +G4double y_77 = y_66; + +G4double x_88 = x_77; +G4double y_88 = y_77 - 8.4*cm; + +G4double x_99 = x_88 - 2.0*cm; +G4double y_99 = y_88; + +G4double x_1010 = x_11 + 8*cos(21*M_PI/180*rad)*cm; +G4double y_1010 = y_11 - 8*sin(21*M_PI/180*rad)*cm; + +std::vector polyd(10); +polyd[0] = G4TwoVector(x_11,y_11); +polyd[1] = G4TwoVector(x_22,y_22); +polyd[2] = G4TwoVector(x_33,y_33); +polyd[3] = G4TwoVector(x_44,y_44); +polyd[4] = G4TwoVector(x_55,y_55); +polyd[5] = G4TwoVector(x_66,y_66); +polyd[6] = G4TwoVector(x_77,y_77); +polyd[7] = G4TwoVector(x_88,y_88); +polyd[8] = G4TwoVector(x_99,y_99); +polyd[9] = G4TwoVector(x_1010,y_1010); + +G4RotationMatrix* rot_l = new G4RotationMatrix(); +rot_l->rotateX(-90*deg); + + +G4ExtrudedSolid* doucone = new G4ExtrudedSolid("DC",polyd,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); +G4LogicalVolume* doucone_log = new G4LogicalVolume(doucone,AlMylar,"DCL",0,0,0,true); +G4VPhysicalVolume* doucone_phy1 = new G4PVPlacement(rot_l,G4ThreeVector(0,pos_s,0),doucone_log,"DCL1",det_log,false,0); +G4VPhysicalVolume* doucone_phy2 = new G4PVPlacement(rot_l,G4ThreeVector(0,-pos_s,0),doucone_log,"DCL2",det_log,false,0); + +G4VisAttributes* dc_c = new G4VisAttributes(); +dc_c->SetColour(1.0,0.0,1.0); +doucone_log->SetVisAttributes(dc_c); + + + +//mirror +G4double xmir1 = x_44; +G4double ymir1 = pos_s; + +G4double xmir2 = x_55; +G4double ymir2 = pos_s; + +G4double xmir3 = x_55; +G4double ymir3 = -pos_s; + +G4double xmir4 = x_44; +G4double ymir4 = -pos_s - thick_ness; + +std::vector polyd1(4); +polyd1[0] = G4TwoVector(xmir1,ymir1); +polyd1[1] = G4TwoVector(xmir2,ymir2); +polyd1[2] = G4TwoVector(xmir3,ymir3); +polyd1[3] = G4TwoVector(xmir4,ymir4); + +G4ExtrudedSolid* doum = new G4ExtrudedSolid("DM",polyd1,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); +G4LogicalVolume* doum_log = new G4LogicalVolume(doum,AlMylar,"DML",0,0,0,true); + +G4RotationMatrix* rot_mir = new G4RotationMatrix(); +rot_mir->rotateY(13.5*deg); + +G4VPhysicalVolume* doum_phy = new G4PVPlacement(rot_mir,G4ThreeVector(thick_ness/2,0,y_11-2*thick_ness),doum_log,"DMP",det_log,false,0); + +G4VisAttributes* mir_r = new G4VisAttributes(); +mir_r->SetColour(0.0,1.0,1.0); +doum_log->SetVisAttributes(mir_r); + + + +//bottom +G4double xbox1 = 4*cm; +G4double ybox1 = (9.35/2)*cm; +G4double zbox1 = thick_ness/2; + +G4Box* bott = new G4Box("Bott",xbox1,ybox1,zbox1); +G4LogicalVolume* bott_log = new G4LogicalVolume(bott,AlMylar,"DCB",0,0,0,true); + +G4RotationMatrix* rot_bott = new G4RotationMatrix(); +rot_bott->rotateY(-21*deg); + +G4VPhysicalVolume* bott_phy = new G4PVPlacement(rot_bott,G4ThreeVector(x_11+ybox1*cos(21*deg),0,-0.5*cm -ybox1*sin(21*deg)),bott_log,"BCP",det_log,false,0); + +bott_log->SetVisAttributes(dc_c); + + +//dconetopbottom +G4double xtop1 = x_66; +G4double ytop1 = pos_s; + +G4double xtop2 = x_77; +G4double ytop2 = pos_s; + +G4double xtop3 = x_77; +G4double ytop3 = -pos_s; + +G4double xtop4 = x_66; +G4double ytop4 = -pos_s; + +std::vector polyd2(4); +polyd2[0] = G4TwoVector(xtop1,ytop1); +polyd2[1] = G4TwoVector(xtop2,ytop2); +polyd2[2] = G4TwoVector(xtop3,ytop3); +polyd2[3] = G4TwoVector(xtop4,ytop4); + +G4ExtrudedSolid* dctop = new G4ExtrudedSolid("DCTB",polyd2,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); +G4LogicalVolume* dctop_log = new G4LogicalVolume(dctop,AlMylar,"DCTBL",0,0,0,true); +G4VPhysicalVolume* dctop_phys1 = new G4PVPlacement(0,G4ThreeVector(0,0,y_55 -thick_ness/2),dctop_log,"DCTBP1",det_log,false,0); +G4VPhysicalVolume* dctop_phys2 = new G4PVPlacement(0,G4ThreeVector(0,0,y_88 +thick_ness/2),dctop_log,"DCTBP2",det_log,false,0); + +dctop_log->SetVisAttributes(dc_c); + + +//cone rear +G4double xrear1 = 1.10*cm; +G4double yrear1 = (9.35/2)*cm; +G4double zrear1 = thick_ness/2; + +G4Box* dcrear = new G4Box("DCREAR",xrear1,yrear1,zrear1); +G4LogicalVolume* dcrear_log = new G4LogicalVolume(dcrear,AlMylar,"DCRL1",0,0,0,true); + +G4RotationMatrix* rot_re = new G4RotationMatrix(); +rot_re->rotateY(90*deg); + +G4VPhysicalVolume* dcrear_phys = new G4PVPlacement(rot_re,G4ThreeVector(11*cm+xbox1*cos(21*deg),0,-7-xbox1*cos(21*deg)),dcrear_log,"DCRP",det_log,false,0); + + + +dcrear_log->SetVisAttributes(dc_c); + + +/* +//Double cut 2cm Moller cone sides +G4double X1 = 3.0*cm; +G4double Y1 = -1.0*cm; + +G4double X2 = -5.0*cm; +G4double Y2 = -1.0*cm; + +G4double X3 = X2; +G4double Y3 = -Y2; + +G4double X4 = X1; +G4double Y4 = Y3; + +G4double X5 = X4 + 8*cos(13.5*M_PI/180*rad)*cm; +G4double Y5 = Y4 + 8*sin(13.5*M_PI/180*rad)*cm; + +G4double X6 = X5 + 2.5*cm; +G4double Y6 = Y5; + +G4double X7 = X6; +G4double Y7 = Y5 - 8.4*cm; + +G4double X8 = X5; +G4double Y8 = Y7; + +G4double X9 = X1 + 8*cos(20.5*M_PI/180*rad)*cm; +G4double Y9 = Y1 - 8*sin(20.5*M_PI/180*rad)*cm; + + + + + + + + +//LGDC +G4double lgsx1 = x_66; //X6; +G4double lgsy1 = y_66; //Y6; + +G4double lgsx2 = x_66 + 2.53*cm; +G4double lgsy2 = y_66; + +G4double lgsx3 = lgsx2 + 40.13*cos(6.5*deg)*cm; +G4double lgsy3 = lgsy2 - 40.13*sin(6.5*deg)*cm; + +G4double lgsx4 = x_66 + 2*cm + 39.6*cos(6.5*deg)*cm; +G4double lgsy4 = lgsy1 - 9.35*cm - 39.6*sin(6.5*deg)*cm; + +G4double lgsx5 = x_66 + 2*cm; +G4double lgsy5 = lgsy1 - 9.35*cm; + +G4double lgsx6 = lgsx1; +G4double lgsy6 = lgsy1 - 9.35*cm; + +std::vector polyd4(6); +polyd4[0] = G4TwoVector(lgsx1,lgsy1); +polyd4[1] = G4TwoVector(lgsx2,lgsy2); +polyd4[2] = G4TwoVector(lgsx3,lgsy3); +polyd4[3] = G4TwoVector(lgsx4,lgsy4); +polyd4[4] = G4TwoVector(lgsx5,lgsy5); +polyd4[5] = G4TwoVector(lgsx6,lgsy6); + +G4ExtrudedSolid* lgdc = new G4ExtrudedSolid("LGDC", polyd4,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); +G4LogicalVolume* lgdc_log = new G4LogicalVolume(lgdc, AlMylar, "LGDCL", 0,0,0,true); + +G4RotationMatrix* rot = new G4RotationMatrix(); +rot->rotateX(-90*deg); + +G4VPhysicalVolume* lgdc_phys1 = new G4PVPlacement(rot,G4ThreeVector(0,pos_s+thick_ness,0),lgdc_log, "LGDCP1", det_log,false,0); +G4VPhysicalVolume* lgdc_phys2 = new G4PVPlacement(rot,G4ThreeVector(0,-pos_s-thick_ness,0),lgdc_log, "LGDCP2", det_log,false,0); + +G4VisAttributes* lg = new G4VisAttributes(); +lg->SetColour(0.0,0.0,1.0); +lgdc_log->SetVisAttributes(lg); + + +std::vectorpolyd5(4); +polyd5[0] = G4TwoVector(lgsx2,pos_s+1.5*thick_ness); +polyd5[1] = G4TwoVector(lgsx2, -pos_s-1.5*thick_ness); +polyd5[2] = G4TwoVector(lgsx3+2*thick_ness, -pos_s-1.5*thick_ness); +polyd5[3] = G4TwoVector(lgsx3+2*thick_ness, pos_s+1.5*thick_ness); + +G4ExtrudedSolid* lgdct = new G4ExtrudedSolid("LGDCT",polyd5,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); +G4LogicalVolume* lgdct_log = new G4LogicalVolume(lgdct, AlMylar,"LGDCTL",0,0,0,true); + +G4RotationMatrix* rot_lgbot = new G4RotationMatrix(); +rot_lgbot->rotateY(-6.5*deg); + +G4VPhysicalVolume* lgdct_phys = new G4PVPlacement(rot_lgbot, G4ThreeVector(-thick_ness,0,5*cm -1.5*thick_ness),lgdct_log,"LGDCTP",det_log,false,0); + +lgdct_log->SetVisAttributes(lg); + +std::vectorpolyd6(4); +polyd6[0] = G4TwoVector(lgsx5,pos_s+1.5*thick_ness); +polyd6[1] = G4TwoVector(lgsx5, -pos_s-1.5*thick_ness); +polyd6[2] = G4TwoVector(lgsx4, -pos_s-1.5*thick_ness); +polyd6[3] = G4TwoVector(lgsx4, pos_s+1.5*thick_ness); + +G4ExtrudedSolid* lgdcb = new G4ExtrudedSolid("LGDCB",polyd6,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); +G4LogicalVolume* lgdcb_log = new G4LogicalVolume(lgdcb, AlMylar,"LGDCTL",0,0,0,true); + +G4VPhysicalVolume* lgdcb_phys = new G4PVPlacement(rot_lgbot, G4ThreeVector(thick_ness,0,y_88+2*thick_ness),lgdcb_log,"LGDCBP",det_log,false,0); + +lgdct_log->SetVisAttributes(lg); +lgdcb_log->SetVisAttributes(lg); + + +std::vectorpolyd7(4); +polyd7[0] = G4TwoVector(lgsx1,pos_s+1.5*thick_ness); +polyd7[1] = G4TwoVector(lgsx2, pos_s+1.5*thick_ness); +polyd7[2] = G4TwoVector(lgsx2, -pos_s-1.5*thick_ness); +polyd7[3] = G4TwoVector(lgsx1, -pos_s-1.5*thick_ness); + +G4ExtrudedSolid* lt = new G4ExtrudedSolid("LT",polyd7,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); +G4LogicalVolume* lt_log = new G4LogicalVolume(lt, AlMylar,"LTL",0,0,0,true); + +G4VPhysicalVolume* lt_phys = new G4PVPlacement(0, G4ThreeVector(0,0,y_55+thick_ness),lt_log,"LTP",det_log,false,0); + +std::vectorpolyd8(4); +polyd8[0] = G4TwoVector(lgsx5,pos_s+1.5*thick_ness); +polyd8[1] = G4TwoVector(lgsx5,-pos_s-1.5*thick_ness); +polyd8[2] = G4TwoVector(lgsx6,-pos_s-1.5*thick_ness); +polyd8[3] = G4TwoVector(lgsx6,pos_s+1.5*thick_ness); + +G4ExtrudedSolid* lb = new G4ExtrudedSolid("LB",polyd8,thick_ness/2,G4TwoVector(0,0),1.0,G4TwoVector(0,0),1.0); +G4LogicalVolume* lb_log = new G4LogicalVolume(lb, AlMylar,"LGDCTL",0,0,0,true); + +G4VPhysicalVolume* lb_phys = new G4PVPlacement(0, G4ThreeVector(0,0,y_88-thick_ness),lb_log,"LGDCBP",det_log,false,0); + +lt_log->SetVisAttributes(lg); +lb_log->SetVisAttributes(lg); + + +//Double cut 2cm Moller cone sides +G4double X1 = 3.0*cm; +G4double Y1 = -1.0*cm; + +G4double X2 = -5.0*cm; +G4double Y2 = -1.0*cm; + +G4double X3 = X2; +G4double Y3 = -Y2; + +G4double X4 = X1; +G4double Y4 = Y3; + +G4double X5 = X4 + 8*cos(13.5*M_PI/180*rad)*cm; +G4double Y5 = Y4 + 8*sin(13.5*M_PI/180*rad)*cm; + +G4double X6 = X5 + 2.5*cm; +G4double Y6 = Y5; + +G4double X7 = X6; +G4double Y7 = Y5 - 8.4*cm; + +G4double X8 = X5; +G4double Y8 = Y7; + +G4double X9 = X1 + 8*cos(20.5*M_PI/180*rad)*cm; +G4double Y9 = Y1 - 8*sin(20.5*M_PI/180*rad)*cm; +*/ + + +//--cathode + + G4double anini=0*deg; + G4double anspan=360.0*deg; + G4double cin = 0; + G4double cout = 7*cm; //big enough to collect all the photons + G4double clngth = 3*mm; + + + + + + + G4Tubs* cath = new G4Tubs("CATH",cin,cout,clngth,anini,anspan); + //---cathode logical + G4LogicalVolume* cath_log = new G4LogicalVolume(cath,CATH,"cath_log",0,0,0); + //---make it sensitive + qsimDetector* cathSD = new qsimDetector("cathSD", 2); //---ID is 2 + + SDman->AddNewDetector(cathSD); + cath_log->SetSensitiveDetector(cathSD); + + G4VisAttributes *cathatt = new G4VisAttributes(); + cathatt->SetColour(1.0, 1.0, 0.2); + cath_log->SetVisAttributes(cathatt); + +//Blackening Tape + G4double blkX,blkY,blkZ; + blkX = 1.0*cm; + blkY = 1.0*cm; + blkZ = 0.0001*cm; + + + G4Box* black = new G4Box("black",blkX,blkY,blkZ); + + G4LogicalVolume* blk_log = new G4LogicalVolume(black,Mylar,"blkMylar",0,0,0); + + + G4VisAttributes *blk_color = new G4VisAttributes(); + blk_color->SetColour(1.0,0.0,0.0); + //place in light guide + blk_log->SetVisAttributes(blk_color); + G4RotationMatrix* rotate_black = new G4RotationMatrix; + rotate_black->rotateY(-11.5*M_PI/180*rad); + G4VPhysicalVolume* black_phys = new G4PVPlacement(rotate_black, G4ThreeVector(-(2+thick_ness)+(xxun-xxdeux)-5*cm,0.0*cm,y_7+8*thick_ness+5.0*cm*sin(11.5*M_PI/180*rad)),blk_log,"BLACK",det_log,false,0); - - - - // Surfaces - - - // quartz - G4OpticalSurface* OpQuartzSurface = new G4OpticalSurface("QuartzSurface"); - OpQuartzSurface->SetType(dielectric_dielectric); - OpQuartzSurface->SetFinish(ground); - OpQuartzSurface->SetModel(glisur); - OpQuartzSurface->SetPolish(fQuartzPolish); - - G4LogicalBorderSurface* QuartzSurface = - new G4LogicalBorderSurface("QuartzSurface",quartz_phys,det_phys,OpQuartzSurface); - - // mirrors and cathode - G4OpticalSurface* MOpSurface = new G4OpticalSurface("MirrorOpSurface"); - G4OpticalSurface* CTHOpSurface = new G4OpticalSurface("CathodeOpSurface"); - - MOpSurface -> SetType(dielectric_metal); - MOpSurface -> SetFinish(polishedlumirrorair); - MOpSurface -> SetModel(glisur); - - CTHOpSurface -> SetType(dielectric_metal); - CTHOpSurface -> SetFinish(polishedlumirrorair); - CTHOpSurface -> SetModel(glisur); - - const G4int num = 2; - G4double Ephoton[num] = {2.038*eV, 4.144*eV}; - - G4MaterialPropertiesTable* MOpSurfaceProperty = new G4MaterialPropertiesTable(); - G4MaterialPropertiesTable* COpSurfaceProperty = new G4MaterialPropertiesTable(); - G4MaterialPropertiesTable* TubeSurfaceProperty = new G4MaterialPropertiesTable(); - - MOpSurfaceProperty -> AddProperty("REFLECTIVITY",PhotonEnergy,Reflectivity1,nEntries); - - MOpSurface -> SetMaterialPropertiesTable(MOpSurfaceProperty); - - COpSurfaceProperty -> AddProperty("REFLECTIVITY",PhotonEnergy,Reflectivity2,nEntries); - COpSurfaceProperty -> AddProperty("EFFICIENCY",PhotonEnergy,EfficiencyArray,nEntries); - - - CTHOpSurface -> SetMaterialPropertiesTable(COpSurfaceProperty); - - if (fDetMode == 0) { - G4LogicalSkinSurface* TubeSurface_1 = new - G4LogicalSkinSurface("TubeMirrorOpS_1",tmirror_log,MOpSurface); - G4LogicalSkinSurface* lightguideSurface = new - G4LogicalSkinSurface("lightguideOps",lightguide_log,MOpSurface); - } - - G4LogicalSkinSurface* CathSurface = new - G4LogicalSkinSurface("CathOpS1", cath_log,CTHOpSurface); - - // Generate & Add Material Properties Table attached to the optical surfaces - - //OpticalQuartzSurface + + + //put cathode in the det box + G4RotationMatrix* rotate_cath =new G4RotationMatrix; +// rotate_cath->rotateY(-96.5*deg); + rotate_cath->rotateY(-101.5*M_PI/180*rad); + + //for moller quartz study +// G4VPhysicalVolume* cath_phys = new G4PVPlacement(rotate_cath,G4ThreeVector(lgsx3 -thick_ness, thick_ness,-4.2*cm -6*thick_ness),cath_log,"Cathode", det_log,false,0); + G4VPhysicalVolume* cath_phys = new G4PVPlacement(rotate_cath,G4ThreeVector(lx3-8*thick_ness,thick_ness,-q_lz - 4*thick_ness),cath_log,"Cathode",det_log,false,0); +// G4VPhysicalVolume* cath_phys = new G4PVPlacement(rotate_cath,G4ThreeVector(lx3-8*thick_ness,thick_ness,-q_lz- 4*thick_ness),cath_log,"Cathode",det_log,false,0); + +//for super elastic quartz study + // G4VPhysicalVolume* cath_phys = new G4PVPlacement(rotate_cath,G4ThreeVector(2.5*cm,0,0),cath_log,"Cathode", det_log,false,0); + + // last step....put det_log into world + G4RotationMatrix* rotate_det =new G4RotationMatrix; + + + G4int step = 1; + G4double worldAngle = /*11.5*deg + */beam_angle; + if(worldAngle < 0) + step = -1; + rotate_det->rotateY(worldAngle); //no rotation corresponding to world system + G4VPhysicalVolume* det_phys = new G4PVPlacement(rotate_det,G4ThreeVector((-(x_6-x_5+thick_ness)+step*(10.305)*sin(((worldAngle/deg)/*-11.5*/)*M_PI/180*rad))*cm,0,0),det_log,"detector_phys",world_log,false,0); + +// world is buit already + + // Surfaces + // quartz optical surface + G4OpticalSurface* OpQuartzSurface = new G4OpticalSurface("QuartzSurface"); + OpQuartzSurface->SetType(dielectric_dielectric); + OpQuartzSurface->SetFinish(ground); + OpQuartzSurface->SetModel(glisur); + OpQuartzSurface->SetPolish(fQuartzPolish); + + const G4int num=2; + G4double Ephoton[num] = {2.038*eV, 4.144*eV}; + + //OpticalQuartzSurface G4double RefractiveIndex[num] = {1.46, 1.46}; G4double SpecularLobe[num] = {0.3, 0.3}; G4double SpecularSpike[num] = {0.2, 0.2}; G4double Backscatter[num] = {0.2, 0.2}; - - G4MaterialPropertiesTable* myST1 = new G4MaterialPropertiesTable(); - - myST1->AddProperty("RINDEX", Ephoton, RefractiveIndex, num); - myST1->AddProperty("SPECULARLOBECONSTANT", Ephoton, SpecularLobe, num); - myST1->AddProperty("SPECULARSPIKECONSTANT", Ephoton, SpecularSpike, num); - myST1->AddProperty("BACKSCATTERCONSTANT", Ephoton, Backscatter, num); - - OpQuartzSurface->SetMaterialPropertiesTable(myST1); - + + G4MaterialPropertiesTable* table_surface_quartz = new G4MaterialPropertiesTable(); + + table_surface_quartz->AddProperty("RINDEX", Ephoton, RefractiveIndex, num); + table_surface_quartz->AddProperty("SPECULARLOBECONSTANT", Ephoton, SpecularLobe, num); + table_surface_quartz->AddProperty("SPECULARSPIKECONSTANT", Ephoton, SpecularSpike, num); + table_surface_quartz->AddProperty("BACKSCATTERCONSTANT", Ephoton, Backscatter, num); + + OpQuartzSurface->SetMaterialPropertiesTable(table_surface_quartz); + + G4LogicalBorderSurface* QuartzSurface = new G4LogicalBorderSurface("QuartzSurface",quartz_phys,det_phys,OpQuartzSurface); + + //reflector surface + G4OpticalSurface* OpReflectorSurface = new G4OpticalSurface("ReflectorSurface"); + OpReflectorSurface->SetType(dielectric_metal); + OpReflectorSurface->SetFinish(polishedlumirrorair); + OpReflectorSurface->SetModel(glisur); + + + G4MaterialPropertiesTable* table_surface_mirror = new G4MaterialPropertiesTable(); + table_surface_mirror->AddProperty("REFLECTIVITY",PhotonEnergy,Reflectivity3,nEntries); + + OpReflectorSurface->SetMaterialPropertiesTable(table_surface_mirror); + + G4LogicalSkinSurface* MirrorSurface = new G4LogicalSkinSurface("ReflectorSurface",ref_log,OpReflectorSurface); + //G4LogicalSkinSurface* MirrorSurface = new G4LogicalSkinSurface("ReflectorSurface",doum_log,OpReflectorSurface); + + G4OpticalSurface* OpMylarSurface = new G4OpticalSurface("MylarSurface"); + OpMylarSurface->SetType(dielectric_metal); + OpMylarSurface->SetFinish(polishedlumirrorair); + OpMylarSurface->SetModel(glisur); + + G4MaterialPropertiesTable* table_surface_mylar = new G4MaterialPropertiesTable(); + table_surface_mylar->AddProperty("REFLECTIVITY",PhotonEnergy,Reflect_LG,nEntries); + + OpMylarSurface->SetMaterialPropertiesTable(table_surface_mylar); + + + G4LogicalSkinSurface* LGSurface1 = new G4LogicalSkinSurface("LGSurface1",l_g_slog,OpMylarSurface); + G4LogicalSkinSurface* LGSurface2 = new G4LogicalSkinSurface("LGSurface2",lgt_log,OpMylarSurface); + G4LogicalSkinSurface* LGSurface3 = new G4LogicalSkinSurface("LGSurface3",lgt_log1,OpMylarSurface); + G4LogicalSkinSurface* LGSurface4 = new G4LogicalSkinSurface("LGSurface4",lgb_log,OpMylarSurface); + G4LogicalSkinSurface* LGSurface5 = new G4LogicalSkinSurface("LGSurface5",lgb_log2,OpMylarSurface); + + G4LogicalSkinSurface* Cone_Side = new G4LogicalSkinSurface("ConeSide",cone_log,OpMylarSurface); + G4LogicalSkinSurface* Cone_bott = new G4LogicalSkinSurface("ConeBottom",coneb_log,OpMylarSurface); + + //G4LogicalSkinSurface* LGSurface1 = new G4LogicalSkinSurface("LGSurface1",lgdct_log,OpMylarSurface); + //G4LogicalSkinSurface* LGSurface2 = new G4LogicalSkinSurface("LGSurface2",lgdcb_log,OpMylarSurface); + //G4LogicalSkinSurface* LGSurface3 = new G4LogicalSkinSurface("LGSurface3",lb_log,OpMylarSurface); + //G4LogicalSkinSurface* LGSurface4 = new G4LogicalSkinSurface("LGSurface4",lt_log,OpMylarSurface); + //4LogicalSkinSurface* LGSurface5 = new G4LogicalSkinSurface("LGSurface5",lgdc_log,OpMylarSurface); + + //G4LogicalSkinSurface* Cone_Side = new G4LogicalSkinSurface("ConeSide",doucone_log,OpMylarSurface); + // G4LogicalSkinSurface* Cone_bott = new G4LogicalSkinSurface("ConeBottom",bott_log,OpMylarSurface); + // G4LogicalSkinSurface* Cone_rear = new G4LogicalSkinSurface("ConeRear",dcrear_log,OpReflectorSurface); + // G4LogicalSkinSurface* Cone_top = new G4LogicalSkinSurface("ConeTop",dctop_log,OpMylarSurface); + + + // cathode surface + + G4OpticalSurface* OpCathodeSurface = new G4OpticalSurface("CathodeSurface"); + OpCathodeSurface->SetType(dielectric_metal); + OpCathodeSurface->SetFinish(polishedlumirrorair); + OpCathodeSurface->SetModel(glisur); + + G4MaterialPropertiesTable* table_surface_cathode = new G4MaterialPropertiesTable(); + table_surface_cathode->AddProperty("REFLECTIVITY", PhotonEnergy, Reflectivity4, nEntries); + table_surface_cathode->AddProperty("EFFICIENCY", PhotonEnergy, Efficiency4, nEntries ); + + OpCathodeSurface->SetMaterialPropertiesTable(table_surface_cathode); + + G4LogicalSkinSurface *cathodeSurface=new G4LogicalSkinSurface("CathodeOpsurface",cath_log,OpCathodeSurface); + //G4LogicalBorderSurface is a class for surfaces defined by the boundary of two physical volumes +//G4LogicalSkinSurface is a class for the surface surronding a single logical volume + +//Blackening Surface +// + G4OpticalSurface* blkSurface = new G4OpticalSurface("blackening"); + blkSurface -> SetType(dielectric_metal); + blkSurface -> SetFinish(polishedlumirrorair); + blkSurface -> SetModel(glisur); + + G4MaterialPropertiesTable* tab_blkSurface = new G4MaterialPropertiesTable(); + tab_blkSurface -> AddProperty("REFLECTIVITY",PhotonEnergy,Reflect_Blk,nEntries); + + blkSurface -> SetMaterialPropertiesTable(tab_blkSurface); + + //always return the physical World return world_phys; -} -//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... - - - - +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/src/qsimMessenger.cc b/src/qsimMessenger.cc index b39736c..bb5bd9e 100644 --- a/src/qsimMessenger.cc +++ b/src/qsimMessenger.cc @@ -11,6 +11,7 @@ #include "qsimEventAction.hh" #include "qsimPrimaryGeneratorAction.hh" #include "qsimSteppingAction.hh" +#include "qsimOpticalPhysics.hh" #include "G4UImanager.hh" #include "G4UIdirectory.hh" @@ -28,6 +29,7 @@ qsimMessenger::qsimMessenger(){ fevact = NULL; fprigen = NULL; fStepAct = NULL; + foptical = NULL; fRemollDir = new G4UIdirectory("/qsim/"); fRemollDir->SetGuidance("UI commands of this code"); @@ -41,6 +43,13 @@ qsimMessenger::qsimMessenger(){ seedCmd->SetParameterName("seed", false); // new + fAllowCerenkovCmd = new G4UIcmdWithABool("/qsim/fAllowCerenkov",this); + fAllowCerenkovCmd->SetGuidance("Turn on Cerenkov"); + fAllowCerenkovCmd->SetParameterName("cerenkov", false); + + fAllowScintillationCmd = new G4UIcmdWithABool("/qsim/fAllowScintillation",this); + fAllowScintillationCmd->SetGuidance("Turn on Scintillation"); + fAllowScintillationCmd->SetParameterName("scintillation", false); //fStandModeCmd = new G4UIcmdWithAnInteger("/qsim/fStandMode",this); //fStandModeCmd->SetGuidance("Set fStandMode to an option"); @@ -58,6 +67,10 @@ qsimMessenger::qsimMessenger(){ fSourceModeCmd->SetGuidance("Set fSourceMode to an option"); fSourceModeCmd->SetParameterName("sourcemode", false); + fLGReflectivityCmd = new G4UIcmdWithADouble("/qsim/fLGReflectivity",this); + fLGReflectivityCmd->SetGuidance("Set fLGReflectivity to a value"); + fLGReflectivityCmd->SetParameterName("lgreflectivity",false); + fQuartzPolishCmd = new G4UIcmdWithADouble("/qsim/fQuartzPolish",this); fQuartzPolishCmd->SetGuidance("Set fQuartzPolish to a value"); fQuartzPolishCmd->SetParameterName("quartzpolish",false); @@ -137,6 +150,14 @@ void qsimMessenger::SetNewValue(G4UIcommand* cmd, G4String newValue){ CLHEP::HepRandom::setTheSeed(seed); } + if (cmd == fAllowCerenkovCmd ) { + foptical->AllowCerenkovSet(fAllowCerenkovCmd->GetNewBoolValue(newValue)); + } + + if (cmd == fAllowScintillationCmd ) { + foptical->AllowScintillationSet(fAllowScintillationCmd->GetNewBoolValue(newValue)); + } + /* if (cmd == fStandModeCmd ) { @@ -159,6 +180,11 @@ void qsimMessenger::SetNewValue(G4UIcommand* cmd, G4String newValue){ fprigen->SourceModeSet(x); } + if (cmd == fLGReflectivityCmd ) { + G4double x = fLGReflectivityCmd->GetNewDoubleValue(newValue); + fdetcon->LGReflectivitySet(x); + } + if (cmd == fQuartzPolishCmd ) { G4double x = fQuartzPolishCmd->GetNewDoubleValue(newValue); fdetcon->fQuartzPolish = x; diff --git a/src/qsimOpticalPhysics.cc b/src/qsimOpticalPhysics.cc index 20d699a..fbecb6a 100644 --- a/src/qsimOpticalPhysics.cc +++ b/src/qsimOpticalPhysics.cc @@ -3,6 +3,22 @@ #include "qsimOpticalPhysics.hh" +qsimOpticalPhysics::qsimOpticalPhysics() +{ + theWLSProcess = NULL; + theScintProcess = NULL; + theCerenkovProcess = NULL; + theBoundaryProcess = NULL; + theAbsorptionProcess = NULL; + theRayleighScattering = NULL; + theMieHGScatteringProcess = NULL; + + fAllowCerenkov = true; + fAllowScintillation = true; + + AbsorptionOn = true; +} + qsimOpticalPhysics::qsimOpticalPhysics(G4bool toggle) : G4VPhysicsConstructor("Optical") { @@ -14,6 +30,9 @@ qsimOpticalPhysics::qsimOpticalPhysics(G4bool toggle) theRayleighScattering = NULL; theMieHGScatteringProcess = NULL; + fAllowCerenkov = true; + fAllowScintillation = true; + AbsorptionOn = toggle; } @@ -30,6 +49,12 @@ void qsimOpticalPhysics::ConstructParticle() } #include "G4ProcessManager.hh" +void qsimOpticalPhysics::AllowCerenkovSet(G4bool set = true) { + fAllowCerenkov = set; +} +void qsimOpticalPhysics::AllowScintillationSet(G4bool set = true) { + fAllowScintillation = set; +} void qsimOpticalPhysics::ConstructProcess() { @@ -100,11 +125,11 @@ FIXME: Add to verbosity responsiveness FatalException,o.str().c_str()); } - if(theCerenkovProcess->IsApplicable(*particle)){ + if(fAllowCerenkov && theCerenkovProcess->IsApplicable(*particle)){ pManager->AddProcess(theCerenkovProcess); pManager->SetProcessOrdering(theCerenkovProcess,idxPostStep); } - if(theScintProcess->IsApplicable(*particle)){ + if(fAllowScintillation && theScintProcess->IsApplicable(*particle)){ pManager->AddProcess(theScintProcess); pManager->SetProcessOrderingToLast(theScintProcess,idxAtRest); pManager->SetProcessOrderingToLast(theScintProcess,idxPostStep); diff --git a/src/qsimPrimaryGeneratorAction.cc b/src/qsimPrimaryGeneratorAction.cc index e3eaa3d..b35384f 100644 --- a/src/qsimPrimaryGeneratorAction.cc +++ b/src/qsimPrimaryGeneratorAction.cc @@ -39,7 +39,7 @@ bool qsimPrimaryGeneratorAction::Thetaspectrum(double Th) { //} // allow user modifications of private member and functional modifiable definition of primary generator variables -void qsimPrimaryGeneratorAction::SourceModeSet(G4int mode = 0) { +void qsimPrimaryGeneratorAction::SourceModeSet(G4int mode = 1) { fSourceMode = mode; // 0 is cosmic mode // 1 is beam mode @@ -67,8 +67,10 @@ void qsimPrimaryGeneratorAction::SourceModeSet(G4int mode = 0) { fEmin = 855.0*MeV; // = 250 MeV at Mainz fEmax = 855.0*MeV; // = 1.063 Gev for JLab - fthetaMin = 0.0*deg; - fthetaMax = 0.0*deg; + //fthetaMin = 0.*deg; + //fthetaMax = 0.*deg; + fthetaMin = 0.*deg; //Trying to make the lightguide perp. to the z axis + fthetaMax = /*9*/0.*deg; } else if (fSourceMode==2){ @@ -89,6 +91,7 @@ qsimPrimaryGeneratorAction::qsimPrimaryGeneratorAction() { fParticleGun = new G4ParticleGun(n_particle); fDefaultEvent = new qsimEvent(); + //fZ = -0.11*m; fZ = -0.52*m; } @@ -158,18 +161,30 @@ void qsimPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) { double randTheta, randPhi; double tanth, tanph; - if (fSourceMode == 0 || fSourceMode == 1) { + if (fSourceMode == 0 /* || fSourceMode == 1*/) { + bool goodTheta = false; + while ( goodTheta == false ) { + randTheta = CLHEP::RandFlat::shoot( fthetaMin, fthetaMax ); + goodTheta = Thetaspectrum(randTheta); + } + + randPhi = CLHEP::RandFlat::shoot( 0.0,360.0)*deg ; + + pX = sin(randTheta)*cos(randPhi)*p; + pY = sin(randTheta)*sin(randPhi)*p; + pZ = cos(randTheta)*p; + } + + if (fSourceMode == 1) { bool goodTheta = false; while ( goodTheta == false ) { randTheta = CLHEP::RandFlat::shoot( fthetaMin, fthetaMax ); goodTheta = Thetaspectrum(randTheta); } - - randPhi = CLHEP::RandFlat::shoot( 0.0,360.0)*deg ; - pX = sin(randTheta)*cos(randPhi)*p; - pY = sin(randTheta)*sin(randPhi)*p; - pZ = cos(randTheta)*p; + pX = sin(randTheta)*p; + pY = 0; + pZ = cos(randTheta)*p; } if (fSourceMode == 2) {