From 7b489612f64d7fc81bc440d3219d1645c990cef8 Mon Sep 17 00:00:00 2001 From: James Szalkie Date: Fri, 29 May 2026 16:48:53 -0400 Subject: [PATCH] more plotting troubleshooting --- Armory/anasenMS.cpp | 12 +- ELoss/E_vs_x_alpha.dat | 1000 ++++++++--------- ELoss/PCEnergyAnalysis.py | 685 +++++++---- .../PCEnergyAnalysis.cpython-310.pyc | Bin 0 -> 28871 bytes ELoss/anasenMS(ElossVersion).cpp | 518 +++++++++ 5 files changed, 1509 insertions(+), 706 deletions(-) create mode 100644 ELoss/__pycache__/PCEnergyAnalysis.cpython-310.pyc create mode 100644 ELoss/anasenMS(ElossVersion).cpp diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 4e42223..9671eba 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -100,13 +100,13 @@ int main(int argc, char **argv){ // then sequential decay of 21Na* -> p + 20Ne TransferReaction transfer; - transfer.SetA(27, 13, 0); // 18Ne projectile - transfer.SetIncidentEnergyAngle((72/27.0), 0, 0); // KEA in MeV/u, theta and phi in degree + transfer.SetA(18, 10, 0); // 18Ne projectile + transfer.SetIncidentEnergyAngle((4.4), 0, 0); // KEA in MeV/u, theta and phi in degree transfer.Seta(4, 2); // 4He target - transfer.Setb(4, 2); // outgoing proton from the primary transfer - transfer.SetB(27, 13); // 21Na* heavy product + transfer.Setb(1, 1); // outgoing proton from the primary transfer + transfer.SetB(21, 11); // 21Na* heavy product - bool enableSequentialDecay = false; // turning to false to disable sequential decay for now, can be set to true to enable + bool enableSequentialDecay = true; // turning to false to disable sequential decay for now, can be set to true to enable const int decayDaughterA = 20; const int decayDaughterZ = 10; const int decayEjectA = 1; @@ -239,7 +239,7 @@ int main(int argc, char **argv){ // reconstructed SX3 hit position double sx3X, sx3Y, sx3Z; - tree1->Branch("sx3X", &sx3X, "sx3X/D"); // reconstructed X position from SX3 (with optional smearing) + tree1->Branch("sx3X", &sx3X, "sx3X/D"); // reconstructed X position from SX3 (with optional smearing) in mm tree1->Branch("sx3Y", &sx3Y, "sx3Y/D"); // reconstructed Y position from SX3 (with optional smearing) tree1->Branch("sx3Z", &sx3Z, "sx3Z/D"); // reconstructed Z position from SX3 (with optional smearing) diff --git a/ELoss/E_vs_x_alpha.dat b/ELoss/E_vs_x_alpha.dat index 3783cea..169f0c3 100644 --- a/ELoss/E_vs_x_alpha.dat +++ b/ELoss/E_vs_x_alpha.dat @@ -1,501 +1,501 @@ Distance_cm Energy_MeV --0.0 40.0 -3.3168136638228574 39.91985971943888 -6.627897879905463 39.83971943887776 -9.933251562889753 39.75957915831663 -13.232873624423497 39.679438877755516 -16.526762973147783 39.59929859719439 -19.8149185146818 39.519158316633266 -23.09733915161019 39.43901803607215 -26.374023783470314 39.358877755511024 -29.644971306736817 39.2787374749499 -32.910180614808766 39.19859719438878 -36.169650597996686 39.11845691382766 -39.42338014350692 39.03831663326653 -42.67136813542852 38.958176352705415 -45.9136134547201 38.87803607214429 -49.15011497919392 38.797895791583166 -52.38087158350261 38.71775551102205 -55.60588213912573 38.63761523046092 -58.8251455143537 38.5574749498998 -62.03866057427424 38.47733466933868 -65.2464261807587 38.39719438877756 -68.44844119244583 38.31705410821643 -71.64470446472787 38.236913827655314 -74.83521484973677 38.15677354709419 -78.01997119632735 38.07663326653307 -81.19897235006462 37.99649298597195 -84.37221715320692 37.91635270541082 -87.53970444469165 37.836212424849705 -90.701433060121 37.75607214428858 -93.85740183174507 37.675931863727456 -97.00760958844737 37.59579158316634 -100.15205515573032 37.515651302605214 -103.29073735569816 37.43551102204409 -106.42365500704213 37.35537074148297 -109.55080692502582 37.27523046092185 -112.67219192146771 37.19509018036072 -115.78780880472625 37.114949899799605 -118.89765637968482 37.03480961923848 -122.00173344773411 36.954669338677355 -125.10003880675694 36.87452905811624 -128.19257125111295 36.79438877755511 -131.27932957162074 36.71424849699399 -134.36031255554238 36.63410821643287 -137.43551898656793 36.553967935871746 -140.50494764479728 36.47382765531062 -143.56859730672446 36.393687374749504 -146.62646674522176 36.31354709418838 -149.67855472952147 36.233406813627255 -152.7248600251998 36.15326653306614 -155.76538139416087 36.07312625250501 -158.80011759461806 35.99298597194389 -161.82906738107778 35.91284569138277 -164.8522295043231 35.832705410821646 -167.86960271139492 35.75256513026052 -170.88118574557538 35.672424849699404 -173.88697734637125 35.59228456913828 -176.88697624949492 35.512144288577154 -179.8811811868474 35.43200400801604 -182.86959088650167 35.35186372745491 -185.852204072683 35.27172344689379 -188.82901946575214 35.19158316633267 -191.80003578218796 35.111442885771545 -194.76525173456793 35.03130260521042 -197.72466603155064 34.9511623246493 -200.67827737785848 34.87102204408818 -203.62608447425757 34.790881763527054 -206.56808601754025 34.710741482965936 -209.50428070050725 34.63060120240481 -212.43466721194756 34.55046092184369 -215.3592442366204 34.47032064128257 -218.2780104552374 34.390180360721445 -221.19096454444184 34.31004008016032 -224.0981051767907 34.2298997995992 -226.99943102073627 34.14975951903808 -229.89494074060525 34.06961923847695 -232.78463299658046 33.989478957915836 -235.66850644468212 33.90933867735471 -238.54655973674684 33.829198396793586 -241.41879152040886 33.74905811623247 -244.28520043908105 33.668917835671344 -247.14578513193376 33.58877755511022 -250.00054423387556 33.5086372745491 -252.84947637553412 33.42849699398798 -255.69258018323467 33.34835671342685 -258.5298542789805 33.268216432865735 -261.3612972804335 33.18807615230461 -264.1869078008924 33.107935871743486 -267.0066844492728 33.02779559118237 -269.82062583008775 32.94765531062124 -272.6287305434252 32.86751503006012 -275.4309971849283 32.787374749499 -278.2274243457751 32.70723446893788 -281.0180106126562 32.62709418837675 -283.80275456775445 32.546953907815634 -286.58165478872445 32.46681362725451 -289.35470984867 32.386673346693385 -292.12191831612313 32.30653306613227 -294.88327875502387 32.22639278557114 -297.63878972469655 32.146252505010025 -300.3884497798303 32.0661122244489 -303.13225747045544 31.98597194388778 -305.8702113419229 31.90583166332666 -308.6023099348818 31.825691382765537 -311.3285517852574 31.745551102204413 -314.0489354242286 31.66541082164329 -316.7634593782064 31.58527054108217 -319.4721221688111 31.505130260521046 -322.17492231284956 31.424989979959925 -324.8718583222932 31.344849699398804 -327.56292870425506 31.26470941883768 -330.2481319609667 31.184569138276558 -332.92746658975597 31.104428857715437 -335.6009310830236 31.024288577154312 -338.26852392822 30.94414829659319 -340.93024360782243 30.86400801603207 -343.5860885993116 30.783867735470945 -346.23605737514794 30.703727454909824 -348.8801484027486 30.623587174348703 -351.5183601444636 30.543446893787582 -354.1506910575522 30.463306613226457 -356.77713959415865 30.383166332665336 -359.3977042012891 30.303026052104215 -362.0123833207872 30.22288577154309 -364.62117538930937 30.14274549098197 -367.22407883830175 30.06260521042085 -369.82109209397515 29.982464929859724 -372.41221357728045 29.902324649298603 -374.99744170388465 29.82218436873748 -377.57677488414606 29.742044088176357 -380.1502115230892 29.661903807615236 -382.7177500203804 29.581763527054115 -385.2793887703029 29.50162324649299 -387.8351261617312 29.42148296593187 -390.3849605781067 29.341342685370748 -392.92889039741226 29.261202404809623 -395.4669139921462 29.181062124248502 -397.99902972929783 29.10092184368738 -400.5252359703214 29.020781563126256 -403.0455310711102 28.940641282565135 -405.55991338197157 28.860501002004014 -408.06838124760054 28.78036072144289 -410.5709330070536 28.70022044088177 -413.06756699372346 28.620080160320647 -415.55828153531246 28.539939879759523 -418.043074953806 28.4597995991984 -420.5219455654469 28.37965931863728 -422.9948916807087 28.299519038076156 -425.46191160426855 28.219378757515035 -427.92300363498157 28.139238476953913 -430.3781660658535 28.05909819639279 -432.8273971840138 27.978957915831668 -435.27069527068915 27.898817635270547 -437.70805860117616 27.818677354709425 -440.13948544481434 27.7385370741483 -442.5649740649586 27.65839679358718 -444.9845227189526 27.57825651302606 -447.3981296581007 27.498116232464934 -449.80579312764064 27.417975951903813 -452.20751136671623 27.33783567134269 -454.60328260834945 27.257695390781567 -456.99310507941243 27.177555110220446 -459.37697700060016 27.097414829659325 -461.75489658640225 27.0172745490982 -464.12686204507446 26.93713426853708 -466.4928715786113 26.856993987975958 -468.85292338271745 26.776853707414833 -471.20701564677915 26.696713426853712 -473.5551465538364 26.61657314629259 -475.89731428055404 26.536432865731467 -478.2335169971932 26.456292585170345 -480.5637528675828 26.376152304609224 -482.88802004909087 26.2960120240481 -485.20631669259507 26.21587174348698 -487.5186409424546 26.135731462925857 -489.82499093648084 26.055591182364733 -492.12536480590785 25.97545090180361 -494.4197606753639 25.89531062124249 -496.7081766628418 25.815170340681366 -498.9906108796695 25.735030060120245 -501.26706143048074 25.654889779559124 -503.53752641318584 25.574749498998 -505.80200391894147 25.494609218436878 -508.0604920321216 25.414468937875757 -510.3129888302875 25.334328657314632 -512.5594923841576 25.25418837675351 -514.8000007575781 25.17404809619239 -517.0345120074929 25.093907815631265 -519.263024183913 25.013767535070144 -521.485535329887 24.933627254509023 -523.7020434814707 24.853486973947902 -525.9125466676966 24.773346693386777 -528.1170429105435 24.693206412825656 -530.3155302249065 24.613066132264535 -532.5080066185662 24.53292585170341 -534.6944700921579 24.45278557114229 -536.874918639141 24.37264529058117 -539.049350245769 24.292505010020044 -541.2177628910574 24.212364729458923 -543.3801545467539 24.1322244488978 -545.5365231773071 24.052084168336677 -547.6868667398351 23.971943887775556 -549.8311831840947 23.891803607214435 -551.9694704524508 23.81166332665331 -554.1017264798438 23.73152304609219 -556.2279491937597 23.651382765531068 -558.3481365141979 23.571242484969943 -560.4622863536397 23.491102204408822 -562.5703966170173 23.4109619238477 -564.6724652016819 23.330821643286576 -566.7684899973718 23.250681362725455 -568.8584688861811 23.170541082164334 -570.9423997425275 23.09040080160321 -573.0202804331207 23.01026052104209 -575.0921088169301 22.930120240480967 -577.1578827451534 22.849979959919843 -579.2176000611835 22.76983967935872 -581.2712586005771 22.6896993987976 -583.3188561910223 22.609559118236476 -585.3603906523058 22.529418837675355 -587.395859796281 22.449278557114233 -589.4252614268355 22.36913827655311 -591.4485933398578 22.288997995991988 -593.4658533232055 22.208857715430867 -595.477039156672 22.128717434869742 -597.4821486119538 22.04857715430862 -599.4811794526174 21.9684368737475 -601.4741294340668 21.88829659318638 -603.4609963035094 21.808156312625254 -605.4417777999238 21.728016032064133 -607.4164716540256 21.647875751503012 -609.3850755882347 21.567735470941887 -611.347587316641 21.487595190380766 -613.3040045449715 21.407454909819645 -615.2543249705558 21.32731462925852 -617.1985462822927 21.2471743486974 -619.1366661606158 21.167034068136278 -621.0686822774599 21.086893787575153 -622.9945922962256 21.006753507014032 -624.9143938717464 20.92661322645291 -626.8280846502527 20.846472945891787 -628.7356622693376 20.766332665330665 -630.6371243579222 20.686192384769544 -632.5324685362203 20.60605210420842 -634.4216924157031 20.5259118236473 -636.3047935990636 20.445771543086178 -638.1817696801814 20.365631262525053 -640.0526182440865 20.28549098196393 -641.917336866923 20.20535070140281 -643.7759231159138 20.125210420841686 -645.6283745493226 20.045070140280565 -647.4746887164185 19.964929859719444 -649.3148631574383 19.88478957915832 -651.1488954035487 19.804649298597198 -652.9767829768102 19.724509018036077 -654.7985233901377 19.644368737474952 -656.6141141472626 19.56422845691383 -658.4235527426952 19.48408817635271 -660.2268366616846 19.403947895791585 -662.0239633801801 19.323807615230464 -663.8149303647915 19.243667334669343 -665.5997350727489 19.16352705410822 -667.3783749518619 19.083386773547097 -669.1508474404794 19.003246492985976 -670.9171499674477 18.923106212424855 -672.6772799520687 18.84296593186373 -674.4312348040572 18.76282565130261 -676.1790119234986 18.68268537074149 -677.9206087008047 18.602545090180364 -679.6560225166695 18.522404809619243 -681.3852507420253 18.44226452905812 -683.1082907379961 18.362124248496997 -684.825139855852 18.281983967935876 -686.5357954369622 18.201843687374755 -688.2402548127476 18.12170340681363 -689.9385153046316 18.04156312625251 -691.6305742239912 17.961422845691388 -693.3164288721068 17.881282565130263 -694.9960765401108 17.801142284569142 -696.6695145089353 17.72100200400802 -698.3367400492591 17.640861723446896 -699.9977504214529 17.560721442885775 -701.6525428755243 17.480581162324654 -703.3011146510609 17.40044088176353 -704.9434629771719 17.32030060120241 -706.5795850724293 17.240160320641287 -708.2094781448071 17.160020040080163 -709.8331393916186 17.07987975951904 -711.4505659994533 16.99973947895792 -713.0617551441111 16.919599198396796 -714.6667039905352 16.839458917835675 -716.2654096927438 16.759318637274554 -717.8578693937587 16.67917835671343 -719.444080225533 16.599038076152308 -721.024039308876 16.518897795591187 -722.5977437533771 16.438757515030062 -724.1651906573258 16.35861723446894 -725.7263771076302 16.27847695390782 -727.2813001797341 16.198336673346695 -728.8299569375288 16.118196392785574 -730.3723444332649 16.038056112224453 -731.9084597074597 15.957915831663328 -733.4382997888018 15.877775551102205 -734.9618616940534 15.797635270541084 -736.4791424279477 15.717494989979961 -737.9901389830849 15.637354709418839 -739.4948483398224 15.557214428857717 -740.9932674661634 15.477074148296595 -742.4853933176395 15.396933867735472 -743.9712228371906 15.31679358717435 -745.45075295504 15.236653306613228 -746.9239805885646 15.156513026052107 -748.3909026421609 15.076372745490984 -749.8515160071054 14.99623246492986 -751.3058175614104 14.91609218436874 -752.7538041696743 14.835951903807617 -754.1954726829248 14.755811623246494 -755.6308199384583 14.675671342685373 -757.0598427596718 14.59553106212425 -758.4825379558876 14.515390781563127 -759.8989023221725 14.435250501002006 -761.3089326391491 14.355110220440883 -762.7126256727994 14.27496993987976 -764.1099781742607 14.19482965931864 -765.5009868796137 14.114689378757516 -766.885648509661 14.034549098196393 -768.2639597696974 13.954408817635272 -769.6359173492704 13.87426853707415 -771.0015179219305 13.794128256513028 -772.3607581449721 13.713987975951905 -773.7136346591625 13.633847695390783 -775.0601440884598 13.553707414829661 -776.4002830397195 13.473567134268539 -777.7340481023883 13.393426853707416 -779.0614358481847 13.313286573146295 -780.3824428307662 13.233146292585172 -781.6970655853826 13.153006012024049 -783.0053006285136 13.072865731462928 -784.3071444574922 12.992725450901805 -785.60259355011 12.912585170340682 -786.8916443642072 12.83244488977956 -788.1742933372436 12.752304609218438 -789.4505368858518 12.672164328657315 -790.7203714053696 12.592024048096194 -791.9837932693544 12.511883767535071 -793.2407988290729 12.43174348697395 -794.4913844129718 12.351603206412827 -795.7355463261227 12.271462925851704 -796.973280849644 12.191322645290583 -798.2045842400963 12.11118236472946 -799.429452728852 12.031042084168337 -800.6478825214355 11.950901803607216 -801.8598697968359 11.870761523046093 -803.0654107067878 11.79062124248497 -804.2645013750199 11.71048096192385 -805.4571378964712 11.630340681362727 -806.6433163364709 11.550200400801604 -807.8230327298825 11.470060120240483 -808.996283080209 11.38991983967936 -810.1630633586577 11.309779559118237 -811.3233695031636 11.229639278557116 -812.477197417368 11.149498997995993 -813.6245429695516 11.06935871743487 -814.7654019915193 10.989218436873749 -815.8997702774344 10.909078156312626 -817.0276435825998 10.828937875751505 -818.1490176221846 10.748797595190382 -819.2638880698923 10.66865731462926 -820.3722505565687 10.588517034068138 -821.4741006687464 10.508376753507015 -822.5694339471231 10.428236472945892 -823.6582458849706 10.348096192384771 -824.740531926471 10.267955911823648 -825.8162874649771 10.187815631262525 -826.8855078411927 10.107675350701404 -827.9481883412708 10.027535070140281 -829.0043241948232 9.947394789579159 -830.0539105728401 9.867254509018037 -831.0969425855137 9.787114228456915 -832.1334152799617 9.706973947895792 -833.1633236378465 9.62683366733467 -834.1866625728848 9.546693386773548 -835.2034269282424 9.466553106212427 -836.2136114738096 9.386412825651304 -837.2172109033501 9.30627254509018 -838.2142198315192 9.22613226452906 -839.2046327907448 9.145991983967937 -840.1884442279628 9.065851703406814 -841.165648501204 8.985711422845693 -842.1362398760216 8.90557114228457 -843.1002125217543 8.825430861723447 -844.0575605076173 8.745290581162326 -845.0082777986116 8.665150300601203 -845.9523582512454 8.58501002004008 -846.8897956090566 8.50486973947896 -847.8205834979296 8.424729458917836 -848.7447154211957 8.344589178356713 -849.6621847545064 8.264448897795592 -850.5729847404713 8.18430861723447 -851.4771084830486 8.104168336673347 -852.3745489416771 8.024028056112225 -853.2652989251382 7.943887775551103 -854.1493510851362 7.863747494989981 -855.0266979095848 7.783607214428859 -855.8973317155845 7.703466933867736 -856.7612446420812 7.623326653306614 -857.618428642189 7.543186372745492 -858.4688754751639 7.46304609218437 -859.3125766980149 7.382905811623247 -860.149523656735 7.302765531062125 -860.9797074771382 7.222625250501003 -861.803119055286 7.14248496993988 -862.6197490474872 7.062344689378758 -863.4295878598543 6.982204408817636 -864.2326256373995 6.902064128256514 -865.028852252653 6.821923847695391 -865.8182572937881 6.741783567134269 -866.6008300522332 6.661643286573147 -867.3765595097566 6.581503006012024 -868.1454343250053 6.501362725450902 -868.9074428194826 6.42122244488978 -869.6625729629483 6.3410821643286575 -870.4108123582267 6.2609418837675355 -871.1521482254104 6.1808016032064135 -871.8865673854444 6.1006613226452915 -872.6140562430837 6.020521042084169 -873.3346007692144 5.940380761523047 -874.0481864825335 5.860240480961925 -874.7547984305863 5.780100200400802 -875.4544211701624 5.69995991983968 -876.1470387470584 5.619819639278558 -876.832634675221 5.539679358717435 -877.5111919152888 5.459539078156313 -878.1826928525622 5.379398797595191 -878.8471192744391 5.299258517034069 -879.5044523473654 5.219118236472946 -880.1546725933622 5.138977955911824 -880.79775986621 5.058837675350702 -881.4336933273827 4.978697394789579 -882.0624514218533 4.898557114228457 -882.6840118539104 4.818416833667335 -883.2983515631593 4.738276553106213 -883.9054467009124 4.65813627254509 -884.5052726072154 4.577995991983968 -885.0978037887983 4.497855711422846 -885.6830138982988 4.4177154308617235 -886.2608757151638 4.3375751503006015 -886.8313611287117 4.2574348697394795 -887.3944411239196 4.177294589178357 -887.9500857706033 4.097154308617235 -888.4982642167707 4.017014028056113 -889.0389446870676 3.93687374749499 -889.5720944873926 3.8567334669338678 -890.0976800169467 3.7765931863727458 -890.6156667891962 3.6964529058116233 -891.1260194634918 3.6163126252505013 -891.6287018893812 3.536172344689379 -892.1236771660092 3.456032064128257 -892.6109077194183 3.3758917835671345 -893.0903554010513 3.295751503006012 -893.5619816113408 3.21561122244489 -894.0257474529553 3.1354709418837676 -894.4816139190865 3.0553306613226456 -894.9295421231245 2.975190380761523 -895.3694935772153 2.8950501002004008 -895.8014305285548 2.8149098196392788 -896.225316363901 2.7347695390781563 -896.6411160947243 2.6546292585170344 -897.0487969377401 2.574488977955912 -897.4483290083564 2.4943486973947895 -897.8396861479281 2.4142084168336675 -898.2228469097687 2.334068136272545 -898.5977957337949 2.253927855711423 -898.9645243456974 2.1737875751503006 -899.3230334239131 2.093647294589178 -899.6733345868216 2.013507014028056 -900.015452764022 1.933366733466934 -900.3494290299925 1.8532264529058118 -900.6753239969484 1.7730861723446896 -900.9932218877797 1.6929458917835674 -901.3032354417734 1.6128056112224451 -901.6055118486687 1.532665330661323 -901.9002399653972 1.4525250501002005 -902.1876591521553 1.3723847695390783 -902.4680701819481 1.292244488977956 -902.7418488488331 1.2121042084168339 -903.0094631542546 1.1319639278557116 -903.2714953361993 1.0518236472945892 -903.5286706029773 0.971683366733467 -903.781895381209 0.8915430861723448 -904.032309435413 0.8114028056112226 -904.2813588341609 0.7312625250501003 -904.5309013671109 0.651122244488978 -904.7833647033009 0.5709819639278558 -905.0419952299258 0.4908416833667335 -905.3112754590231 0.4107014028056113 -905.5976925055561 0.330561122244489 -905.9113738277266 0.25042084168336676 -906.2705145196951 0.17028056112224452 -906.7180616804087 0.09014028056112225 -907.5636766533732 0.01 +-0.0 22.0 +1.0979679668260327 21.955931863727454 +2.1940971205166475 21.911863727454907 +3.2883870865884766 21.867795591182364 +4.380837489806898 21.823727454909818 +5.471447954183834 21.779659318637272 +6.5602181029763385 21.73559118236473 +7.647147558685182 21.691523046092183 +8.732235943052554 21.64745490981964 +9.81548287706108 21.603386773547093 +10.896887980931519 21.559318637274547 +11.976450874121337 21.515250501002004 +13.05417117532327 21.471182364729458 +14.130048502463108 21.42711422845691 +15.204082472698252 21.38304609218437 +16.276272702416268 21.338977955911822 +17.346618807232662 21.294909819639276 +18.41512040198943 21.250841683366733 +19.481777100753597 21.206773547094187 +20.546588516814978 21.16270541082164 +21.60955426268472 21.118637274549098 +22.670673950093835 21.07456913827655 +23.72994718999086 21.03050100200401 +24.787373592540813 20.986432865731462 +25.842952767122856 20.942364729458916 +26.896684322328806 20.898296593186373 +27.948567865961643 20.854228456913827 +28.99860300503325 20.81016032064128 +30.046789345762907 20.766092184368738 +31.09312649357579 20.72202404809619 +32.1376140531007 20.677955911823645 +33.18025162816852 20.633887775551102 +34.22103882181073 20.589819639278556 +35.25997523625709 20.54575150300601 +36.29706047293413 20.501683366733467 +37.33229413246358 20.45761523046092 +38.36567581466003 20.413547094188377 +39.39720511852974 20.36947895791583 +40.4268816422683 20.325410821643285 +41.45470498325903 20.281342685370742 +42.480674738071414 20.237274549098196 +43.50479050245879 20.19320641282565 +44.52705187135674 20.149138276553106 +45.54745843888153 20.10507014028056 +46.56600979832772 20.061002004008014 +47.58270554216661 20.01693386773547 +48.59754526204458 19.972865731462925 +49.61052854878075 19.928797595190378 +50.62165499236534 19.884729458917835 +51.63092418195805 19.84066132264529 +52.63833570588556 19.796593186372746 +53.6438891516403 19.7525250501002 +54.647584105877975 19.708456913827654 +55.649420154415886 19.66438877755511 +56.64939688223124 19.620320641282564 +57.64751387345871 19.576252505010018 +58.64377071138879 19.532184368737475 +59.638166978466 19.48811623246493 +60.630702256286504 19.444048096192383 +61.62137612559633 19.39997995991984 +62.61018816628963 19.355911823647293 +63.59713795740624 19.311843687374747 +64.58222507712985 19.267775551102204 +65.56544910278625 19.223707414829658 +66.54680961084074 19.179639278557115 +67.52630617689674 19.13557114228457 +68.50393837569314 19.091503006012022 +69.47970578110247 19.04743486973948 +70.45360796612907 19.003366733466933 +71.42564450290652 18.959298597194387 +72.39581496269575 18.915230460921844 +73.36411891588315 18.871162324649298 +74.33055593197797 18.82709418837675 +75.29512557961039 18.78302605210421 +76.25782742652956 18.738957915831662 +77.21866103960099 18.694889779559116 +78.17762598480454 18.650821643286573 +79.13472182723244 18.606753507014027 +80.08994813108653 18.562685370741484 +81.04330445967659 18.518617234468937 +81.99479037541754 18.47454909819639 +82.94440543982736 18.43048096192385 +83.89214921352504 18.386412825651302 +84.83802125622775 18.342344689378756 +85.7820211267487 18.298276553106213 +86.72414838299507 18.254208416833666 +87.66440258196505 18.21014028056112 +88.60278327974575 18.166072144288577 +89.53929003151092 18.12200400801603 +90.47392239151802 18.077935871743485 +91.40667991310599 18.033867735470942 +92.33756214869297 17.989799599198395 +93.26656864977318 17.945731462925853 +94.19369896691498 17.901663326653306 +95.11895264975777 17.85759519038076 +96.04232924700958 17.813527054108217 +96.9638283064446 17.76945891783567 +97.8834493749001 17.725390781563124 +98.80119199827392 17.68132264529058 +99.71705572152199 17.637254509018035 +100.63104008865504 17.59318637274549 +101.54314464273607 17.549118236472946 +102.45336892587771 17.5050501002004 +103.36171247923892 17.460981963927853 +104.26817484302227 17.41691382765531 +105.17275555647123 17.372845691382764 +106.0754541578667 17.32877755511022 +106.97627018452452 17.284709418837675 +107.87520317279196 17.24064128256513 +108.77225265804475 17.196573146292586 +109.66741817468426 17.15250501002004 +110.56069925613379 17.108436873747493 +111.4520954348356 17.06436873747495 +112.34160624224785 17.020300601202404 +113.22923120884083 16.976232464929858 +114.11496986409384 16.932164328657315 +114.99882173649193 16.88809619238477 +115.88078635352208 16.844028056112222 +116.76086324166984 16.79995991983968 +117.63905192641596 16.755891783567133 +118.51535193223233 16.71182364729459 +119.38976278257883 16.667755511022044 +120.26228399989908 16.623687374749498 +121.1329151056169 16.579619238476955 +122.00165562013264 16.53555110220441 +122.86850506281881 16.491482965931862 +123.73346295201638 16.44741482965932 +124.5965288050309 16.403346693386773 +125.45770213812796 16.359278557114227 +126.31698246652924 16.315210420841684 +127.17436930440847 16.271142284569137 +128.02986216488668 16.22707414829659 +128.88346056002803 16.183006012024048 +129.7351640008355 16.138937875751502 +130.58497199724593 16.09486973947896 +131.43288405812586 16.050801603206413 +132.27889969126642 16.006733466933866 +133.12301840337872 15.96266533066132 +133.96523970008897 15.918597194388775 +134.8055630859335 15.87452905811623 +135.64398806435372 15.830460921843684 +136.4805141376908 15.78639278557114 +137.31514080718068 15.742324649298595 +138.14786757294857 15.69825651302605 +138.97869393400356 15.654188376753504 +139.80761938823292 15.61012024048096 +140.63464343239667 15.566052104208415 +141.4597655621217 15.521983967935869 +142.28298527189583 15.477915831663324 +143.10430205506194 15.43384769539078 +143.9237154038118 15.389779559118235 +144.7412248091799 15.345711422845689 +145.55682976103697 15.301643286573144 +146.37052974808364 15.2575751503006 +147.18232425784382 15.213507014028053 +147.99221277665777 15.169438877755509 +148.80019478967552 15.125370741482964 +149.60626978084966 15.08130260521042 +150.41043723292825 15.037234468937873 +151.2126966274474 14.993166332665329 +152.013047444724 14.949098196392784 +152.81148916384805 14.905030060120238 +153.60802126267467 14.860961923847693 +154.40264321781652 14.816893787575149 +155.19535450463545 14.772825651302604 +155.98615459723433 14.728757515030058 +156.77504296844853 14.684689378757513 +157.56201908983738 14.640621242484968 +158.34708243167537 14.596553106212422 +159.13023246294298 14.552484969939878 +159.9114686513177 14.508416833667333 +160.69079046316452 14.464348697394788 +161.46819736352637 14.420280561122242 +162.2436888161142 14.376212424849697 +163.01726428329718 14.332144288577153 +163.7889232260923 14.288076152304606 +164.55866510415385 14.244008016032062 +165.32648937576295 14.199939879759517 +166.0923954978164 14.155871743486973 +166.8563829258156 14.111803607214426 +167.61845111385506 14.067735470941882 +168.3785995146108 14.023667334669337 +169.1368275793284 13.97959919839679 +169.89313475781069 13.935531062124246 +170.64752049840533 13.891462925851702 +171.39998424799208 13.847394789579157 +172.1505254519697 13.80332665330661 +172.89914355424253 13.759258517034066 +173.64583799720694 13.715190380761522 +174.39060822173732 13.671122244488975 +175.13345366717175 13.62705410821643 +175.8743737712974 13.582985971943886 +176.61336797033562 13.538917835671342 +177.3504356989266 13.494849699398795 +178.08557639011363 13.45078156312625 +178.81878947532732 13.406713426853706 +179.550074384369 13.36264529058116 +180.279430545394 13.318577154308615 +181.00685738489457 13.27450901803607 +181.73235432768223 13.230440881763526 +182.45592079686992 13.18637274549098 +183.17755621385342 13.142304609218435 +183.89725999829275 13.09823647294589 +184.61503156809283 13.054168336673344 +185.33087033938364 13.0101002004008 +186.04477572650026 12.966032064128255 +186.75674714196214 12.92196392785571 +187.46678399645188 12.877895791583164 +188.17488569879364 12.83382765531062 +188.88105165593106 12.789759519038075 +189.58528127290452 12.745691382765528 +190.28757395282784 12.701623246492984 +190.98792909686472 12.65755511022044 +191.68634610420426 12.613486973947895 +192.3828243720361 12.569418837675348 +193.0773632955248 12.525350701402804 +193.7699622677839 12.48128256513026 +194.46062067984911 12.437214428857713 +195.1493379206508 12.393146292585168 +195.83611337698616 12.349078156312624 +196.5209464334904 12.30501002004008 +197.20383647260732 12.260941883767533 +197.88478287455916 12.216873747494988 +198.56378501731592 12.172805611222444 +199.24084227656365 12.128737474949897 +199.9159540256721 12.084669338677353 +200.5891196356617 12.040601202404808 +201.26033847516945 11.996533066132264 +201.92960991041446 11.952464929859717 +202.59693330516203 11.908396793587173 +203.26230802068744 11.864328657314628 +203.92573341573856 11.820260521042082 +204.58720884649745 11.776192384769537 +205.24673366654142 11.732124248496993 +205.90430722680273 11.688056112224448 +206.55992887552756 11.643987975951902 +207.21359795823383 11.599919839679357 +207.86531381766818 11.555851703406812 +208.51507579376167 11.511783567134266 +209.16288322358452 11.467715430861722 +209.80873544129977 11.423647294589177 +210.4526317781158 11.379579158316632 +211.09457156223758 11.335511022044086 +211.73455411881673 11.291442885771541 +212.37257876990066 11.247374749498997 +213.00864483437996 11.20330661322645 +213.64275162793493 11.159238476953906 +214.27489846298053 11.115170340681361 +214.90508464861014 11.071102204408817 +215.53330949053787 11.02703406813627 +216.15957229103935 10.982965931863726 +216.78387234889124 10.938897795591181 +217.40620895930917 10.894829659318635 +218.02658141388406 10.85076152304609 +218.64498900051703 10.806693386773546 +219.2614310033526 10.762625250501001 +219.8759067027103 10.718557114228455 +220.4884153750143 10.67448897795591 +221.09895629272197 10.630420841683366 +221.70752872424987 10.58635270541082 +222.3141319338984 10.542284569138275 +222.9187651817746 10.49821643286573 +223.52142772371266 10.454148296593186 +224.12211881119285 10.41008016032064 +224.72083769125817 10.366012024048095 +225.3175836064291 10.32194388777555 +225.91235579461627 10.277875751503004 +226.50515348903053 10.23380761523046 +227.09597591809145 10.189739478957915 +227.684822305333 10.14567134268537 +228.27169186930715 10.101603206412824 +228.85658382348484 10.057535070140279 +229.43949737615483 10.013466933867734 +230.02043173031976 9.969398797595188 +230.59938608358956 9.925330661322644 +231.17635962807248 9.881262525050099 +231.7513515502631 9.837194388777554 +232.32436103092783 9.793126252505008 +232.89538724498715 9.749058116232463 +233.4644293613954 9.704989979959919 +234.0314865430171 9.660921843687373 +234.5965579465005 9.616853707414828 +235.15964272214777 9.572785571142283 +235.72074001378203 9.528717434869739 +236.2798489586111 9.484649298597192 +236.83696868708753 9.440581162324648 +237.3920983227656 9.396513026052103 +237.94523698215434 9.352444889779557 +238.496383774567 9.308376753507012 +239.04553780196673 9.264308617234468 +239.5926981588085 9.220240480961923 +240.13786393187684 9.176172344689377 +240.6810342001195 9.132104208416832 +241.22220803447732 9.088036072144288 +241.7613844977093 9.043967935871741 +242.2985626442135 8.999899799599197 +242.83374151984367 8.955831663326652 +243.36692016172083 8.911763527054108 +243.89809759804044 8.867695390781561 +244.4272728478745 8.823627254509017 +244.95444492096888 8.779559118236472 +245.4796128175355 8.735490981963926 +246.00277552803905 8.691422845691381 +246.52393203297873 8.647354709418837 +247.04308130266418 8.603286573146292 +247.56022229698598 8.559218436873746 +248.07535396518023 8.515150300601201 +248.5884752455873 8.471082164328656 +249.0995850654046 8.42701402805611 +249.60868234043286 8.382945891783566 +250.1157659748164 8.338877755511021 +250.6208348607766 8.294809619238476 +251.12388787833876 8.25074148296593 +251.62492389505215 8.206673346693385 +252.12394176570305 8.16260521042084 +252.62094033202058 8.118537074148295 +253.115918422375 8.07446893787575 +253.60887485146878 8.030400801603205 +254.0998084200197 7.98633266533066 +254.58871791443607 7.942264529058115 +255.07560210648393 7.89819639278557 +255.56045975294577 7.854128256513025 +256.0432895952711 7.81006012024048 +256.52409035921784 7.765991983967934 +257.002860754485 7.72192384769539 +257.47959947433674 7.677855711422844 +257.95430519521597 7.6337875751503 +258.42697657634955 7.589719438877754 +258.8976122593431 7.54565130260521 +259.3662108677658 7.501583166332664 +259.8327710067253 7.457515030060119 +260.2972912624318 7.413446893787574 +260.7597702017518 7.369378757515029 +261.22020637175046 7.325310621242484 +261.6785982992233 7.281242484969939 +262.1349444902161 7.237174348697394 +262.5892434295332 7.193106212424849 +263.04149358023403 7.149038076152303 +263.49169338311737 7.1049699398797586 +263.9398412561931 7.060901803607213 +264.3859355941415 7.0168336673346685 +264.82997476775915 6.972765531062123 +265.2719571233923 6.9286973947895785 +265.71188098235626 6.884629258517033 +266.14974464034105 6.8405611222444875 +266.5855463668035 6.796492985971943 +267.01928440434455 6.7524248496993975 +267.45095696807255 6.708356713426853 +267.8805622449515 6.6642885771543074 +268.3080983931342 6.620220440881763 +268.73356354128055 6.576152304609217 +269.1569557878599 6.532084168336672 +269.57827320043776 6.488016032064127 +269.9975138149468 6.443947895791582 +270.41467563494103 6.399879759519037 +270.82975663083425 6.355811623246492 +271.2427547391214 6.311743486973947 +271.6536678615831 6.267675350701402 +272.06249386447337 6.223607214428856 +272.4692305776899 6.179539078156312 +272.87387579392725 6.135470941883766 +273.2764272678121 6.091402805611222 +273.6768827150216 6.047334669338676 +274.07523981138314 6.003266533066132 +274.4714961919573 5.959198396793586 +274.8656494501023 5.915130260521041 +275.25769713652073 5.871062124248496 +275.64763675828914 5.826993987975951 +276.035465777869 5.782925851703406 +276.4211816121007 5.738857715430861 +276.8047816311802 5.694789579158316 +277.1862631576179 5.650721442885771 +277.565623465181 5.606653306613225 +277.9428597778192 5.562585170340681 +278.3179692685742 5.518517034068135 +278.69094905847334 5.4744488977955905 +279.0617962154081 5.430380761523045 +279.4305077529982 5.3863126252505005 +279.7970806294419 5.342244488977955 +280.1615117463532 5.2981763527054095 +280.52379794758724 5.254108216432865 +280.88393601805495 5.2100400801603195 +281.2419226825278 5.165971943887775 +281.59775460443467 5.1219038076152295 +281.95142838465176 5.077835671342685 +282.30294056028777 5.033767535070139 +282.65228760346605 4.989699398797594 +282.99946592010576 4.945631262525049 +283.34447184870487 4.901563126252504 +283.68730165912706 4.857494989979959 +284.02795155139574 4.813426853707414 +284.3664176544985 4.769358717434869 +284.7026960252052 4.725290581162324 +285.03678264690404 4.681222444889778 +285.36867342845926 4.637154308617234 +285.6983642030958 4.593086172344688 +286.0258507273159 4.549018036072144 +286.35112867985276 4.504949899799598 +286.67419366066827 4.460881763527054 +286.9950411900013 4.416813627254508 +287.3136667074735 4.372745490981963 +287.63006557126135 4.328677354709418 +287.9442330573434 4.284609218436873 +288.2561643588316 4.240541082164328 +288.56585458539826 4.196472945891783 +288.8732987628097 4.152404809619238 +289.1784918325797 4.108336673346693 +289.4814286517556 4.064268537074147 +289.78210399285376 4.020200400801603 +290.08051254395946 3.9761322645290575 +290.37664890900965 3.9320641282565125 +290.67050760827874 3.887995991983967 +290.9620830790882 3.843927855711422 +291.25136967676355 3.799859719438877 +291.5383616758648 3.755791583166332 +291.8230532717176 3.711723446893787 +292.10543858227584 3.667655310621242 +292.3855116503494 3.623587174348697 +292.66326644623234 3.5795190380761515 +292.9386968707727 3.5354509018036064 +293.21179675892574 3.4913827655310614 +293.4825598838393 3.4473146292585164 +293.75097996152226 3.4032464929859714 +294.0170506561529 3.3591783567134263 +294.2807655860894 3.3151102204408813 +294.54211833064954 3.271042084168336 +294.80110243773447 3.226973947895791 +295.0577114323765 3.182905811623246 +295.31193882630066 3.138837675350701 +295.56377812859665 3.0947695390781558 +295.81322285760746 3.0507014028056108 +296.0602665541514 3.0066332665330657 +296.30490279620517 2.9625651302605203 +296.54712521518843 2.9184969939879752 +296.7869275140029 2.87442885771543 +297.02430348699534 2.830360721442885 +297.25924704202885 2.78629258517034 +297.4917522248668 2.742224448897795 +297.72181324609204 2.69815631262525 +297.9494245108078 2.6540881763527047 +298.1745806513904 2.6100200400801596 +298.3972765635919 2.5659519038076146 +298.61750744632036 2.5218837675350696 +298.8352688454589 2.4778156312625246 +299.05055670212374 2.4337474949899796 +299.26336740579967 2.3896793587174345 +299.4736978528418 2.345611222444889 +299.6815455108801 2.301543086172344 +299.88690848972414 2.257474949899799 +300.0897856194288 2.213406813627254 +300.29017653625476 2.169338677354709 +300.48808177734 2.125270541082164 +300.68350288499187 2.081202404809619 +300.87644252161334 2.0371342685370735 +301.06690459639896 1.993066132264529 +301.2548944050708 1.9489979959919836 +301.4404187840838 1.9049298597194386 +301.6234862809111 1.8608617234468936 +301.80410734222795 1.8167935871743486 +301.9822945220607 1.7727254509018033 +302.15806271225074 1.7286573146292583 +302.33142939792094 1.6845891783567133 +302.5024149410324 1.640521042084168 +302.6710428955916 1.596452905811623 +302.83734035863887 1.552384769539078 +303.0013383618331 1.508316633266533 +303.16307230927816 1.4642484969939877 +303.3225824682487 1.4201803607214427 +303.47991452071557 1.3761122244488977 +303.63512018510113 1.3320440881763524 +303.7882579195974 1.2879759519038074 +303.93939372075243 1.2439078156312624 +304.08860203401105 1.1998396793587174 +304.23596679666525 1.1557715430861721 +304.381582638461 1.1117034068136271 +304.52555627124974 1.067635270541082 +304.66800810699596 1.0235671342685368 +304.80907415375975 0.9794989979959918 +304.9489082528044 0.9354308617234468 +305.08768473791486 0.8913627254509017 +305.2256016220649 0.8472945891783566 +305.3628844492392 0.8032264529058115 +305.4997909942559 0.7591583166332665 +305.6366170566244 0.7150901803607214 +305.7737036849217 0.6710220440881762 +305.9114463007515 0.6269539078156312 +306.05030639112084 0.5828857715430861 +306.1908267489405 0.538817635270541 +306.3336517435569 0.4947494989979959 +306.47955495034967 0.45068136272545084 +306.6294779701976 0.40661322645290576 +306.7845870883725 0.3625450901803607 +306.9463600739776 0.3184769539078156 +307.1167276661543 0.27440881763527053 +307.2983234411945 0.23034068136272542 +307.49497395910805 0.18627254509018035 +307.7128080105875 0.14220440881763527 +307.9631635274653 0.09813627254509016 +308.2668481468979 0.054068136272545086 +308.76087549467394 0.01 diff --git a/ELoss/PCEnergyAnalysis.py b/ELoss/PCEnergyAnalysis.py index 28f9b79..9af6182 100644 --- a/ELoss/PCEnergyAnalysis.py +++ b/ELoss/PCEnergyAnalysis.py @@ -11,30 +11,54 @@ from scipy.interpolate import interp1d import uproot import pycatima as catima from scipy.integrate import cumulative_trapezoid -import matplotlib -matplotlib.use("Agg") +#matplotlib.use("Agg") import matplotlib.pyplot as plt import cmd import shlex import textwrap -import readline -import atexit import os import periodictable as pt import re -from matplotlib.colors import PowerNorm +from matplotlib.colors import LinearSegmentedColormap +import matplotlib as mpl -#Run program from terminal or IDE, and prompts will provide user steps +# ROOT-like styling +plt.rcParams.update({ -#histfile = os.path.expanduser("~/.uproot_shell_history") + # Figure + "figure.facecolor": "white", + "axes.facecolor": "white", -#try: -# readline.read_history_file(histfile) -#except FileNotFoundError: -# pass + # ROOT-like axes + "axes.linewidth": 1.2, -#atexit.register(readline.write_history_file, histfile) + # Ticks on all sides + "xtick.top": True, + "ytick.right": True, + # Tick direction inward + "xtick.direction": "in", + "ytick.direction": "in", + + # Major/minor ticks + "xtick.major.size": 10, + "ytick.major.size": 10, + "xtick.minor.size": 5, + "ytick.minor.size": 5, + + # Smaller ROOT-like fonts + "font.size": 12, +}) + +base_cmap = plt.cm.jet + +# Stop before the red region +colors = base_cmap(np.linspace(0.2, 0.65, 256)) + +yellow_jet = LinearSegmentedColormap.from_list( + 'yellow_jet', + colors +) #data = [z, mass_u, maximum MeV, name] alpha_data = [2, 4.0026, 40, "alpha"] proton_data = [1, 1.0078, 20, "proton"] @@ -172,7 +196,7 @@ def get_interpolators(particle): x, E, bounds_error=False, - fill_value=0.0 + fill_value="extrapolate" ) x_of_E = interp1d( @@ -258,43 +282,67 @@ def update_plot_data(name, values): return plot_data.append((name, values)) -def process_file(filename): +def calculate_distance_tree1(vx, vy, vz, sx3x, sx3y, sx3z): + """Calculate 3D distance from vertex to SX3 hit position for tree1""" + dx = (sx3x - vx) / 10 + dy = (sx3y - vy) / 10 + dz = (sx3z - vz) / 10 + return np.sqrt(dx*dx + dy*dy + dz*dz) - tree = uproot.open(filename)["tree"] +def calculate_distance_tree2(vz, theta, z_max=34.86): + """Calculate distance along trajectory from vertex to z_max for tree2 + theta is polar angle from z-axis in radians""" + cos_theta = np.cos(theta) + cos_theta = np.where(np.abs(cos_theta) < 1e-10, 1e-10, cos_theta) + return np.abs((z_max - vz) / cos_theta) - branches = ["Tb", "thetab", "sx3Z"] - data = tree.arrays(branches, library="np") +def load_tree_arrays(tree, treename, max_events=None): + branches = ["Tb", "thetab", "sx3Z", "vX", "vY", "vZ"] + if treename == 'tree1': + branches.extend(["sx3X", "sx3Y"]) + return tree.arrays(branches, library="np", entry_stop=max_events) + + +def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86): + data = load_tree_arrays(tree, treename, max_events) Ei = data["Tb"] theta = np.radians(data["thetab"]) - sx3Z = data["sx3Z"] + vX = data["vX"] + vY = data["vY"] + vZ = data["vZ"] - mask = np.sin(theta) != 0 + if treename == 'tree1': + sx3X = data["sx3X"] + sx3Y = data["sx3Y"] + sx3Z = data["sx3Z"] + mask = ~np.isnan(sx3X) & ~np.isnan(sx3Y) & ~np.isnan(sx3Z) + else: + sx3Z = np.full_like(vZ, z_max) # Assume sx3Z is at z_max for tree2 + mask = ~np.isnan(Ei) & ~np.isnan(theta) & ~np.isnan(vZ) Ei = Ei[mask] theta = theta[mask] - sx3Z = sx3Z[mask] + vX = vX[mask] + vY = vY[mask] + vZ = vZ[mask] + + if treename == 'tree1': + sx3X = sx3X[mask] + sx3Y = sx3Y[mask] + sx3Z = sx3Z[mask] + dsx3 = calculate_distance_tree1(vX, vY, vZ, sx3X, sx3Y, sx3Z) + else: + dsx3 = calculate_distance_tree2(vZ, theta, z_max=z_max) sin_theta = np.sin(theta) - - radii = np.array([3.2, 4.2, 6.6]) - + sin_theta = np.where(sin_theta != 0, sin_theta, 1e-10) + radii = np.array([3.7, 4.2]) dA = radii[0] / sin_theta dC = radii[1] / sin_theta - dsx3 = radii[2] / sin_theta - # Determine particle from filename - lower = filename.lower() - - if "proton" in lower: - particle = "proton" - elif "alpha" in lower: - particle = "alpha" - else: - particle = "proton" - - print(f"Computing energies for {particle}...") + print(f"Computing energies for {particle} ({treename})...") EA = energy_loss(particle, Ei, dA) EC = energy_loss(particle, Ei, dC) @@ -314,6 +362,21 @@ def process_file(filename): "Eprop": Eprop } + +def process_file(filename, treename, particle=None, max_events=None): + tree = uproot.open(filename)[f"{treename}"] + + if particle is None: + lower = filename.lower() + if "proton" in lower: + particle = "proton" + elif "alpha" in lower: + particle = "alpha" + else: + particle = "proton" + print(f"File {filename} particle {particle}, tree {treename}") + return prepare_tree_data(tree, treename, particle, max_events=max_events) + class MyInteractiveApp(cmd.Cmd): def __init__(self): super().__init__() @@ -484,7 +547,7 @@ class MyInteractiveApp(cmd.Cmd): plt.tight_layout() filename = f"Energy_Loss_Curve_{label}.png" plt.savefig(filename, dpi=300, bbox_inches="tight") - plt.close() + plt.show() print(f"Saved plot: {filename}") @@ -601,7 +664,7 @@ class MyInteractiveApp(cmd.Cmd): print("Available trees:", keys) - treeName = arg if len(arg) > 0 else "tree" + treeName = arg if len(arg) > 0 else "tree1" try: # uproot automatically resolves ";1" @@ -696,18 +759,20 @@ class MyInteractiveApp(cmd.Cmd): self.do_set_tree("") print(f"Using TTree: {self.tree}") - branches = [ - "Tb", - "thetab", - "sx3Z" - ] + treename = self.tree.name if hasattr(self.tree, 'name') else "tree1" + + + branches = ["Tb", "thetab", "sx3Z", "vX", "vY", "vZ", "sx3X", "sx3Y"] + Tb_key = "Tb" + theta_key = "thetab" + if max_events: n_events = max_events else: n_events = self.tree.num_entries - print(f"Loading {n_events} events...") + print(f"Loading {n_events} events from {treename}...") data = self.tree.arrays( branches, @@ -715,25 +780,42 @@ class MyInteractiveApp(cmd.Cmd): entry_stop=max_events ) global Ei - Ei = data["Tb"] + Ei = data[Tb_key] global sx3Z sx3Z = data["sx3Z"] - theta = np.radians(data["thetab"]) - - # Remove theta = 0 events - mask = np.sin(theta) != 0 + theta = np.radians(data[theta_key]) + vX = data["vX"] + vY = data["vY"] + vZ = data["vZ"] + + if treename == 'tree1': + sx3X = data["sx3X"] + sx3Y = data["sx3Y"] + mask = ~np.isnan(sx3X) & ~np.isnan(sx3Y) & ~np.isnan(sx3Z) + else: + mask = ~np.isnan(Ei) & ~np.isnan(theta) Ei = Ei[mask] theta = theta[mask] - + sx3Z = sx3Z[mask] + vZ_masked = vZ[mask] + sin_theta = np.sin(theta) + sin_theta = np.where(sin_theta != 0, sin_theta, 1e-10) - radii = np.array([3.2, 4.2, 6.6]) - + if treename == 'tree1': + vX_masked = vX[mask] + vY_masked = vY[mask] + sx3X_masked = sx3X[mask] + sx3Y_masked = sx3Y[mask] + dsx3 = calculate_distance_tree1(vX_masked, vY_masked, vZ_masked, sx3X_masked, sx3Y_masked, sx3Z) + else: + dsx3 = calculate_distance_tree2(vZ_masked, theta, z_max=34.86) + + radii = np.array([3.7, 4.2]) dA = radii[0] / sin_theta dC = radii[1] / sin_theta - dsx3 = radii[2] / sin_theta print("Calculating energy losses...") global EA @@ -830,162 +912,265 @@ class MyInteractiveApp(cmd.Cmd): print("Finished writing augmented ROOT file") def do_make_plots(self, arg): - + import os - + global plot_data - plot_data = [] - + args = shlex.split(arg) - particle = args[0] if len(args) > 0 else "proton" max_events = int(args[1]) if len(args) > 1 else None - + if self.tree is None: self.do_set_tree("") print(f"Using TTree: {self.tree}") - - branches = ["Tb", "thetab", "sx3Z"] - - data = self.tree.arrays( - branches, - library="np", - entry_stop=max_events - ) - - Ei = data["Tb"] - theta = np.radians(data["thetab"]) + + treename = self.tree.name if hasattr(self.tree, 'name') else "tree1" + + data = prepare_tree_data(self.tree, treename, particle, max_events=max_events) + + Ei = data["Ei"] sx3Z = data["sx3Z"] - - mask = np.sin(theta) != 0 - - Ei = Ei[mask] - theta = theta[mask] - sx3Z = sx3Z[mask] - - update_plot_data(f"{particle}_Ei", Ei) - update_plot_data(f"{particle}_sx3Z", sx3Z) - - sin_theta = np.sin(theta) - - radii = np.array([3.2, 4.2, 6.6]) - - dA = radii[0] / sin_theta - dC = radii[1] / sin_theta - dsx3 = radii[2] / sin_theta - - print("Computing energies...") - - EA = energy_loss(particle, Ei, dA) - EC = energy_loss(particle, Ei, dC) - Esx3 = energy_loss(particle, Ei, dsx3) - - Elost = Ei - Esx3 - Eprop = EA - EC - - update_plot_data(f"{particle}_EA", EA) - update_plot_data(f"{particle}_EC", EC) - update_plot_data(f"{particle}_Esx3", Esx3) - update_plot_data(f"{particle}_Elost", Elost) - update_plot_data(f"{particle}_Eprop", Eprop) - - base = f"{particle}_plots" + EA = data["EA"] + EC = data["EC"] + Esx3 = data["Esx3"] + Elost = data["Elost"] + Eprop = data["Eprop"] + + update_plot_data(f"{particle}_{treename}_Ei", Ei) + update_plot_data(f"{particle}_{treename}_sx3Z", sx3Z) + update_plot_data(f"{particle}_{treename}_EA", EA) + update_plot_data(f"{particle}_{treename}_EC", EC) + update_plot_data(f"{particle}_{treename}_Esx3", Esx3) + update_plot_data(f"{particle}_{treename}_Elost", Elost) + update_plot_data(f"{particle}_{treename}_Eprop", Eprop) + + base = f"{particle}_{treename}_plots" os.makedirs(base, exist_ok=True) - - print(f"Saving plots to folder: {base}") - - # 1. Histogram: energy loss + + print(f"Saving plots to folder: {base} ({treename})") + + branch_names = [] + for key in self.tree.keys(): + if isinstance(key, bytes): + branch_names.append(key.decode("ascii")) + elif hasattr(key, "name"): + branch_names.append(key.name) + else: + branch_names.append(str(key)) + branch_names = list(dict.fromkeys(branch_names)) + + if branch_names: + print(f"Creating histograms for {len(branch_names)} branches...") + all_branches = self.tree.arrays(branch_names, library="np", entry_stop=max_events) + + for branch in branch_names: + values = all_branches[branch] + try: + values = np.asarray(values, dtype=float) + except Exception: + print(f"Skipping non-numeric branch: {branch}") + continue + if values.ndim != 1: + print(f"Skipping non-scalar branch: {branch}") + continue + values = values[~np.isnan(values)] + if values.size == 0: + print(f"Skipping empty branch: {branch}") + continue + + plt.figure(figsize=(7,5)) + plt.hist(values, bins=100) + plt.xlabel(branch) + plt.ylabel("Counts") + plt.title(f"{particle} ({treename}) {branch} distribution") + plt.grid(True) + plt.tight_layout() + safe_name = re.sub(r"[^0-9A-Za-z_-]", "_", branch) + plt.savefig(f"{base}/{safe_name}_hist.png", dpi=300) + plt.close() + else: + print("No branches found to histogram.") + plt.figure(figsize=(7,5)) plt.hist(Elost, bins=200) plt.xlabel("Energy Loss (MeV)") plt.ylabel("Counts") - plt.title(f"{particle} Total Energy Loss Distribution") + plt.title(f"{particle} ({treename}) Total Energy Loss Distribution") plt.grid(True) plt.tight_layout() plt.savefig(f"{base}/Elost_hist.png", dpi=300) - plt.close() - - #Histogram: sx3Z + plt.show() + plt.figure(figsize=(7,5)) - plt.hist(sx3Z, bins=100) - plt.xlabel("SX3 Z") + plt.hist(Ei, bins=100, histtype='step', label='Vertex Energy') + plt.hist(Esx3[Esx3 != 0], bins=100, histtype='step', label='SX3 Energy') + plt.xlabel("Energies (MeV)") plt.ylabel("Counts") - plt.title(f"{particle} SX3 Position Distribution") + plt.title(f"{particle} ({treename}) Vertex and SX3 Energy") plt.grid(True) plt.tight_layout() - plt.savefig(f"{base}/sx3Z_hist.png", dpi=300) - plt.close() + plt.savefig(f"{base}/Vertex_vs_SX3_Energy.png", dpi=300) + plt.legend(loc='upper right') + plt.show() + + try: + plt.figure(figsize=(7,5)) + plt.hist(sx3Z, bins=100) + plt.xlabel("SX3 Z") + plt.ylabel("Counts") + plt.title(f"{particle} ({treename}) SX3 Position Distribution") + plt.grid(True) + plt.tight_layout() + plt.savefig(f"{base}/sx3Z_hist.png", dpi=300) + plt.show() + except ValueError: + print("Value error") + + try: + plt.figure(figsize=(7,6)) + plt.hist2d(sx3Z, Elost, bins=200) + plt.ylabel("Energy Loss (MeV)") + plt.xlabel("SX3 Z") + plt.title(f"{particle} ({treename}) Energy Loss vs SX3 Position") + plt.colorbar(label="Counts") + plt.tight_layout() + plt.savefig(f"{base}/Elost_vs_sx3Z.png", dpi=300) + plt.show() + except ValueError: + print("Value error") + + try: + plt.figure(figsize=(7,6)) + plt.hist2d(EA, Esx3, bins=200) + plt.xlabel("EA (MeV)") + plt.ylabel("Esx3 (MeV)") + plt.title(f"{particle} ({treename}) Anode vs SX3 Energy") + plt.colorbar(label="Counts") + plt.tight_layout() + plt.savefig(f"{base}/EA_vs_Esx3.png", dpi=300) + plt.show() + except ValueError: + print("Value error") + + try: + plt.figure(figsize=(7,6)) + plt.hist2d(sx3Z, Eprop, bins=200) + plt.ylabel("EA - EC (MeV)") + plt.xlabel("SX3 Z") + plt.title(f"{particle} ({treename}) Energy Propagation Difference vs Position") + plt.colorbar(label="Counts") + plt.tight_layout() + plt.savefig(f"{base}/Eprop_vs_sx3Z.png", dpi=300) + plt.show() + except ValueError: + print("Value error") + + try: + mask1 = ~np.isnan(Esx3) & (Esx3 > 0) + plt.figure(figsize=(7,6)) + plt.hist2d(Esx3[mask1], Eprop[mask1], bins=200) + plt.ylabel("PCEnergy") + plt.xlabel("SX3 Energy (MeV)") + plt.title(f"{particle} ({treename}) Energy Propagation Difference vs SX3 Energy") + plt.colorbar(label="Counts") + plt.xlim(0,30) + plt.ylim(0,0.45) + plt.tight_layout() + plt.savefig(f"{base}/Eprop_vs_Esx3.png", dpi=300) + plt.show() + except ValueError: + print("Value error") - #2D: Elost vs sx3Z - plt.figure(figsize=(7,6)) - plt.hist2d(sx3Z, Elost, bins=200) - plt.ylabel("Energy Loss (MeV)") - plt.xlabel("SX3 Z") - plt.title(f"{particle} Energy Loss vs SX3 Position") - plt.colorbar(label="Counts") - plt.tight_layout() - plt.savefig(f"{base}/Elost_vs_sx3Z.png", dpi=300) - plt.close() - - #Anode energy vs sx3 energy - plt.figure(figsize=(7,6)) - plt.hist2d(EA, Esx3, bins=200) - plt.xlabel("EA (MeV)") - plt.ylabel("Esx3 (MeV)") - plt.title(f"{particle} Anode vs SX3 Energy") - plt.colorbar(label="Counts") - plt.tight_layout() - plt.savefig(f"{base}/EA_vs_Esx3.png", dpi=300) - plt.close() - - #Prop counter energy loss - plt.figure(figsize=(7,6)) - plt.hist2d(Eprop, sx3Z, bins=200) - plt.xlabel("EA - EC (MeV)") - plt.ylabel("SX3 Z") - plt.title(f"{particle} Energy Propagation Difference vs Position") - plt.colorbar(label="Counts") - plt.tight_layout() - plt.savefig(f"{base}/Eprop_vs_sx3Z.png", dpi=300) - plt.close() - - plt.figure(figsize=(7,6)) - plt.hist2d(Esx3, Eprop, bins=200) - plt.ylabel("PCEnergy") - plt.xlabel("SX3 Energy (MeV)") - plt.title(f"{particle} Energy Propagation Difference vs sx3 Energy") - plt.colorbar(label="Counts") - plt.xlim(0,30) - plt.ylim(0,.45) - plt.tight_layout() - plt.savefig(f"{base}/Eprop_vs_Esx3.png", dpi=300) - plt.close() - - print("Plotting complete.") + print(f"Plotting complete ({treename}).") + def do_hist_comp(self, arg): + particle = arg + self.do_set_tree("tree1") + data1 = prepare_tree_data(self.tree, "tree1", particle, max_events=None) + self.do_set_tree("tree2") + data2 = prepare_tree_data(self.tree, "tree2", particle, max_events=None) + + base = f"{particle}_comp_plots" + os.makedirs(base, exist_ok=True) + + branch_names = [] + for key in self.tree.keys(): + if isinstance(key, bytes): + branch_names.append(key.decode("ascii")) + elif hasattr(key, "name"): + branch_names.append(key.name) + else: + branch_names.append(str(key)) + branch_names = list(dict.fromkeys(branch_names)) + + print(branch_names) + + if branch_names: + print(f"Creating histograms for {len(branch_names)} branches...") + self.do_set_tree("tree1") + all_branches1 = self.tree.arrays(branch_names, library="np", entry_stop=None) + self.do_set_tree("tree2") + all_branches2 = self.tree.arrays(branch_names, library="np", entry_stop=None) + for branch in branch_names: + values1 = all_branches1[branch] + values2 = all_branches2[branch] + try: + values1 = np.asarray(values1, dtype=float) + values2 = np.asarray(values2, dtype=float) + except Exception: + print(f"Skipping non-numeric branch: {branch}") + continue + if values1.ndim != 1 or values2.ndim != 1: + print(f"Skipping non-scalar branch: {branch}") + continue + values1 = values1[~np.isnan(values1)] + values2 = values2[~np.isnan(values2)] + if values1.size == 0 or values2.size == 0: + print(f"Skipping empty branch: {branch}") + continue + + plt.figure(figsize=(7,5)) + plt.hist(values1, bins=100, histtype='step') + plt.hist(values2, bins=100, histtype='step') + plt.xlabel(branch) + plt.ylabel("Counts") + plt.title(f"{particle} {branch} distribution") + plt.grid(True) + plt.tight_layout() + safe_name = re.sub(r"[^0-9A-Za-z_-]", "_", branch) + plt.savefig(f"{base}/{safe_name}_hist.png", dpi=300) + plt.show() + else: + print("No branches found to histogram.") + + def do_dual_plotter(self, arg): args = shlex.split(arg) if len(args) < 2: - print("Usage: make dual plots proton_data.root alpha_data.root") - return + try: + args = ['SimAnasenProton.root', 'SimAnasenAlpha.root'] + except: + print("Usage: make dual plots proton_data.root alpha_data.root") + return file1 = args[0] file2 = args[1] - - #Helper function - - - #Process both files - data1 = process_file(f"../Armory/{file1}") - data2 = process_file(f"../Armory/{file2}") - + #tree1_name = args[2] if len(args) > 2 else 'tree2' + #tree2_name = args[3] if len(args) > 3 else 'tree1' outdir = "dual_plots" + + data1 = process_file(os.path.join("..", "Armory", file1), "tree1") + data2 = process_file(os.path.join("..", "Armory", file2), "tree1") + + #print(f"File one {file1} ({tree1_name}) length {len(data1['Ei'])} \nFile two {file2} ({tree2_name}) length {len(data2['Ei'])}") os.makedirs(outdir, exist_ok=True) print(f"Saving plots to: {outdir}") + #Overlay histogram: Elost plt.figure(figsize=(8,6)) @@ -1015,7 +1200,7 @@ class MyInteractiveApp(cmd.Cmd): plt.tight_layout() plt.savefig(f"{outdir}/Elost_overlay.png", dpi=300) - plt.close() + plt.show() #Overlay histogram: sx3Z plt.figure(figsize=(8,6)) @@ -1046,7 +1231,7 @@ class MyInteractiveApp(cmd.Cmd): plt.tight_layout() plt.savefig(f"{outdir}/sx3Z_overlay.png", dpi=300) - plt.close() + plt.show() #Side-by-side 2D plots fig, axes = plt.subplots(1, 2, figsize=(14,6)) @@ -1077,7 +1262,7 @@ class MyInteractiveApp(cmd.Cmd): plt.tight_layout() plt.savefig(f"{outdir}/Elost_vs_sx3_comparison.png", dpi=300) - plt.close() + plt.show() #EA vs Esx3 overlay scatter plt.figure(figsize=(8,6)) @@ -1106,48 +1291,148 @@ class MyInteractiveApp(cmd.Cmd): plt.tight_layout() plt.savefig(f"{outdir}/EA_vs_Esx3_overlay.png", dpi=300) - plt.close() + plt.show() #PCE vs SiE + mask1 = data1["Esx3"] > 0 + mask2 = data2["Esx3"] > 0 combined_Esx3 = np.concatenate([ - data1["Esx3"], - data2["Esx3"] + (data1["Esx3"][mask1]), + data2["Esx3"][mask2] ]) - combined_Eprop = np.concatenate([ - data1["Eprop"], - data2["Eprop"] + data1["Eprop"][mask1], + data2["Eprop"][mask2] ]) - plt.figure(figsize=(7,6)) - + plt.figure(figsize=(8,6)) plt.hist2d( combined_Esx3, combined_Eprop, - bins=100, - range=[[0,30],[0,0.45]], - cmap='viridis' # same default matplotlib style + bins=200 ) - - plt.ylabel("PCEnergy") plt.xlabel("SX3 Energy (MeV)") - - plt.title( - f'{data1["particle"]} + {data2["particle"]} ' - 'Energy Propagation Difference vs SX3 Energy' - ) - + plt.ylabel("PCEnergy") + plt.xlim(0,35) + plt.ylim(0,0.6) plt.colorbar(label="Counts") - - plt.xlim(0,30) - plt.ylim(0,0.45) - plt.tight_layout() + plt.savefig(f"{outdir}/Eprop_vs_Esx3.png", dpi=300) + plt.show() + + + + plt.figure(figsize=(12,6), facecolor='white') + plt.hist2d( + combined_Esx3, + combined_Eprop, + bins=[500,500], + #range=[[0,35],[0,0.6]], + cmap=yellow_jet, + cmin=1 + ) + plt.xlabel("SX3 Energy (MeV)") + plt.ylabel("PCEnergy") + plt.xlim(0,35) + plt.ylim(0,0.6) + cbar = plt.colorbar() + cbar.set_label("Counts") + plt.minorticks_on() + plt.grid(False) + plt.tight_layout() + plt.savefig( + f"{outdir}/ROOT_style_plot.png", + dpi=300, + facecolor='white' + ) + plt.show() + if False: + plt.figure(figsize=(7,6)) + plt.hist2d( + combined_Esx3, + combined_Eprop, + bins=100, + range=[[0,30],[0,0.45]], + cmap='viridis' + ) + + def power_fit_and_plot(x, y, label, color=None): + + x = np.array(x) + y = np.array(y) + + mask = ( + np.isfinite(x) & + np.isfinite(y) & + (x > 0) & + (y > 0) + ) + + x = x[mask] + y = y[mask] + + logx = np.log(x) + logy = np.log(y) + + b, loga = np.polyfit(logx, logy, 1) + + a = np.exp(loga) + + y_pred = a * x**b + + ss_res = np.sum((y - y_pred)**2) + ss_tot = np.sum((y - np.mean(y))**2) + + r2 = 1 - (ss_res / ss_tot) + + x_fit = np.linspace(np.min(x), np.max(x), 1000) + y_fit = a * x_fit**b + + plt.plot( + x_fit, + y_fit, + linewidth=3, + color=color, + label=f'{label} fit ($R^2$={r2:.4f})' + ) + + print(f"\n{label} power-law fit:") + print(f"y = {a:.6f} * x^{b:.6f}") + print(f"R^2 = {r2:.6f}") + + return a, b, r2 + + + power_fit_and_plot( + data1["Esx3"]["Esx3" > 0], + data1["Eprop"]["Esx3" > 0], + label=data1["particle"], + color='red' + ) + + power_fit_and_plot( + data2["Esx3"], + data2["Eprop"], + label=data2["particle"], + color='cyan' + ) + + plt.ylabel("PCEnergy") + plt.xlabel("SX3 Energy (MeV)") + plt.title( + f'{data1["particle"]} + {data2["particle"]} ' + 'Energy Propagation Difference vs SX3 Energy' + ) + plt.colorbar(label="Counts") + #plt.xlim(0,30) + #plt.ylim(0,0.45) + plt.legend() + plt.tight_layout() + plt.savefig(f"{outdir}/Combined_Eprop_vs_Esx3_fit.png",dpi=300) + plt.show() - plt.savefig(f"{outdir}/Combined_Eprop_vs_Esx3.png", dpi=300) - plt.close() print("Dual plotting complete.") diff --git a/ELoss/__pycache__/PCEnergyAnalysis.cpython-310.pyc b/ELoss/__pycache__/PCEnergyAnalysis.cpython-310.pyc new file mode 100644 index 0000000000000000000000000000000000000000..3380a2dbf26af73a46b97606a11659b539118148 GIT binary patch literal 28871 zcmbV#378z&U0-)~^?gjwp*eJxZuPD-T1k7@+Fc#ij_t+XS$Q?yW7g~5R#(laM?KxW zRW*{PsIeWiYXh=w&@mhg!nS>zQ>->KI zSKZUyGtw^Z{`$RF@2YqF-~auu_qv(SXCnB#_W6m0PhO8izQBj>&j21y;pcoh8i|;Z z*@zi6wR&_`6F)Yq!;jVDi;3ApRHY?nlaiL2O~cpgnZ@jEHcZPc=4T6Gcwlxw@)c)` z;t$Raia#_vB>wR1i1?dkH^GnBM;FIt$HKDXvlCWfVe{%1Gl6f?OyQfxH-m2$-<+9$ zC3hnljUd<7g>AFjt?ib!y2C7(17`6?>}qt@u*Oz*nuC^cLz~?-A2EkM5uM!)e;EEA z_#^Q5nww_#nWMA&&9T`7)^2N$wbvYfB{o~CMP~0ZC(O++MrQ9ex0qYuA2hd_+uAN(Wce)9nQqh`sx3;r?l0rQ}F=*8&lz2;lYN%Js5 z_nB`s?=g=cbiaArJZc_8=z(`c%zMrIUW}OcSpy4ixf+{&t99JG|DtC9w)NJlk<|xp zL@+kmIr)_y^d8Ebv8{@0nMSi=JO%Hu%Cd3v9^=Td}&WR+dJ(z3i}jqzjet~IQR zea@O+w2<38)2uh`#Y!uPU%pUtt@&f}`O3+eK&v%w654LY>Ez&)a$a$+Hos(BlXI1- zRb^4HP`P3`>z?a)%jKHsUhvXa+*gIMPrc7%IbJA?Qre2*s^_qHh@l@EDou^7Mo9osd zqQfeyNjUA#Rs=W$XYg}eF6-z@($rSsW^5(Fupc+|k;qEQPx(ndF%sDlSxKAmm5iBK z$(qTPoSE`7BM~$0XU&YC8;OiX_*=>QX+MwnEbm zlo~M8@N;I-9Gr{I$5sa1=&I%q+=#9e{o-oOFL3E1PWPjJmAM(fG5Ba0wkD4R?D0(vP=KOJg$lPR(UW=_v_!HU4%7mL= zE%+0N9eYJXJ2so+e(d`r=EU{b$`*g~!hpZU+)Tu9ErxzX%x(Ub?~h#9WE42$a^!z0;>3D+x5J7^5!@2319fvCv8F|)_1bjsJXRDy`%-%M{15+ zX;iIpbuq|HH7t96x%`;*WRRRHUvkP<%D41Mv*O;0rr(a-F~S*!%9$AG15)E^|nHV2urc5{JixNZenqZOpE_MH~+Ly++R`A)fB zyJ!VztFg3b*(e=EX9C@{%qu~%V$VBG+qGk8d!V(^ zFGws_9H+b#q^!ltmGYu>DTvoA=Pk6N(Q2Z7PLMcTo_YLslbfZ#Vc52RKYdOd+RxVm)HB-#xydRFP=}~@Dv0PMz$F-a` zq^0y!G!@mqo=Ie)Ijs=O#EM!$8;Yf3Be9f*w3w!8QGOF&PbBmhL#fwualI2;NTf9T zJ$OZxfYkoz2oMQT>S780)B%8QDI@Vu zSniVDlc`*>?aJ~&V@d@imlu1VbN8+2?rA~Hz|XtTfgo9JE^_M&5*MtBX>p=1)-bP^ z9kqmX`yOUXVLvtrfSA5|vl^rTj)36Sa!C(ToV*5sRLAuBX*R$YgwVG&LHEOO_YOAXU0TejV_gAB%_UcOYRFIn~$WG}^noFYZa)k^h(CEc^R zt(Fo&s#UQ;2GlM48APm&U67b6H|NS&C0FQA^^fJ?yoKc?1P)u#w4I+pd1ntCLU<1A zpbg`I@~>&hTr{U3{%g4;hPaB1RU(@9hgg;HEkn#X=K#ns##VF;1BP@Z?#F!%n`Od} z`3XEHW$OdQd;^=ptJ?Hy?aE$_$qVcfB7!84RBdr-v6Ql(=d(6d3pBHCSMXq;mzXPv zxiW{8x&D4;6(%Xyn~w9}koPO3!Eze4Wr5VXHyRH5-m=F9}We1?FtfX~OAgHL*%(^zXvMhNG|h*oIdUGU^FG zJ`O?w;~t#`$;A0o!yT)fcS=baYSlJ2KZurPs8y4t@qJAm2%Bt%yX0S?ltZ9UgI&;U z2M;}~l8t5#0hvESz>lm(ZbVk2_<_X&ZGjj!hIKmV#?(4>_0_nk30;I$YES!7&_=Nv z(bXi8Bq$<2>o*8&F>l;WFJwULkV+v4lF|VPF68jbn*c?ucz`0Plh~a86tI9EaSO|j zMM1=Zav3m_#0jKazAbt=;thtiMffh%nD7&*DIqmQy2>wh%BQ-^m(tTg)bU0SJX<;B zojx@C+})*T553QM?&N{@nRk~2r>>DQk;sWtx1xu;V=-3=rh2h7x3LBdvkx}DbvG}{W|*p90gg801MTxtcGYNb_kD|K)LL5#$S z)^b4^)GeZQPy){K;`wGhNB~`*N7cIQ(dnrxRjWltBgi~SOf^N+RZ0dr5xM=Xn9QIM z)=>^?h@K9TVHD^p)uBLVOh5P+ZOdua!D_YJeGBzDL@Gnrmw}zIXBY6*v>f;8m=4U8 z1D47q5cdsDD@KR#%&`4Y6si(ewm&5PJ1gaFcyK{%cgScBXrK5B;0Eq;%#Zj{pz^4% zJwy%@Ao7APt4HEVW;4#_j)FCG5MGdO)thda8(EO3v|3ifEXC~)GbJG?-Lc=nhrD$K z3#HbWFUx@0h+8)bOD$A`L1J2*f)vgs_8y8v|}JqZdgZ#K3(qo261^`kVdYu#Jt%7J-UYj zO7P0W#<8sIPa=5xQJh371W3_OU<`ebwOB+0erz>D>S!K#7@ROk-Waj3nYb1c$_QkP zHgg^<3j2RvCVv(n@iczV;5Um-aB3ytX>p&CPQR)J$$IU)U9p#g43L$*Ty}sZK~N!h0m91ZAC#gZI8jE?!Op=pbP)>`qQ$8F zlN=5LQTsE209T-FYM(=X?!;i_V^?F`e9R~|A+-ZzBb(qxLXiWvrW8dF@YC?)@G}^| z1P1UFo|Aaa!cW1+_GPBw=YfeJP`D8tptt}c1%^(ka9)L13XO53>3#;n11&6(=#)ry zmB_S901FgQdjfcXeiC?qehLB(`f1cfzumeI`&oI;`Ve^_KIiAf&%1GyWNxP9uO?Oo z8c8>S=cKQj8Hf`o_yDfK2ow=Y`9(8(F$On>g~ph)nKugrnw3F+(9D9|E&_uMAv_R< zhuw@>lx=*_k76?)x*n&P!Nr!3UJPNxN5b$h!kfbI0K%h49|Z+B=8yQB;Ezj*JfA26 z@yE^KS5g4Zk!kw}F`4!c(fKr;SLu9)&JWZ15jy{h&ZpqK8nr*cFs40X{{Wm&DwRKF zk20zV#~VA-Tm;_=+;0fYRI?mG1h+KffTx!xCnwv4;{}#sy`XGQ<98;|r+^=W_!P)l z7aka`Rw*YysL((%#19ZXIE_jpNY5>Sm!QZ%@Sv<-$R=tHdzID3WqsIt_#RTir3w4T z8U8Gt=jr?e9qH37L(-xf4D}60?>GpRa1s0*Tbst~6w(ReY#yn1^u|16ljX2~md@ws z{2ZO1r_)DU4z_FyqyR;)Sc=?Rk8vSBjjSL5v;r2Re94;QyASNk9>j;SozO_Mpk=fn zEL?64qagq3GoL0bd>1GA|Lk8!+%|wU4g#@w^J)cg#n(2>K#8)*2v$K1vbh$_EP;5jK!7!H(_EkDT6Or zO-X(}rGW$$)MU)0NLeGRS*9XBCD4>H)1fcMPPEwRVZ0DBA$2nM#yw%Se* z@dnP-79m7)tj3W^LX)@s&>cVC#G}(QQ%{^eGxO+^Q^vV>PQClx#<`i(Ps~g`;!Vsf zw=84-1*_iLFX)bfqY9)Rlip^P1d7<*Pa+aKw*iowhv2f51241GqFkvAaYfB7mxJW_ zrRqh?b?ifEY7jjeN1PV`vcd7NXl2oP zH456os?Whsm&?^UWu#?=iq3LSeh1D9p`WUVOf;b-0Qd0m%L6+SE+qfw;HtoftsNHWn?Yt|01bk`!;%q2H=5**3@- z143%o@FHH__VzUiu>@iH=8*Di_v1ip2JmzVsJT1e^ zGH-}Uy)4t-rqHTN={zu{oh{c@1~n*&pFz>8mz@%g4HOVv50aAW%4T?>WYg%CY0lh`;d5_?f}ypPcKZaw0BciR8t7Dx1v>S9eJGpyM8T4)#5o6%Iirn$P(Ba^5!~Z$ zmzZvMNb|_M%K%EWe zRZVpxh7R<%Uq{zxein`Y3g1xByiMzygdDeRkNv+<3@tb)JCmV=mu$iXfP#<~ZK#Hv z1>}aPJv*)YT@<$el1}eV7-Z0_a2|gFjkEs-dAuzfH7(5bh7H`(H$-QR1K@TX2WBCq zPJ3i$-l#nsrnA4)X~~9zbhf7rWMMk;39@s-#!%>`QQ1ZL1cm1ixPuUtur6UpC=gGv z4P_M6;6TmgC+2l*-?~fkGy!5o#OV;6>j*=hfel_rP%}-E)BDaTlEvg(LXt#wZkV>L zu`(cN&RfPj=;j8A8|cZG7b}Ld?10AWC^4t7BC#vVx@(ht#zWo-;~@&CjU@?)&d8Wj zz=xPlM#_IqC}NlVxJcb^azphB)7?$k?xR}eR3m8x&s506BX$CPJ zR3owF?TSD{4%weTnjz&fSO@RGnF|nCLXZI&Ue|J;G~Vk}=ApJAPIQf9_uS)+&KMV9 zqH!K_I1pA8a0?ad?K)exz-fT^Cx}#7NzH7pTS}2W!(N`HT$a>yPj8KT+1az-WM|K= zbyn8C;LO{0_Rr9~fJlIPlG}H%7ZG#=LA!%nH|$`0``*`cV152-D+{{9j0-0+RNM%-Yl;z&gB%C)v_e-c? zg6@|P%-a0aYVn4sOvAWA6C41Dfh_Aab1;X?^DPo8T#v70QSZ==$igs;;s*To(+EX7 z?j{O$DAc(YgYj8X-Z3gdk1;=m-DMo68acFa45oCjMm4vHxmjj)Ttpw5xz*1p!qk}zHsJJNicyf$cD_qAd&(=bNtET!h$|+`1iS4zU!R z`_f5R#iY(RuQ={9L_Fz^0j=16Lf?~y>aF2m-&llw5tw{7G#I3!66?E;y=A^3W=r~c zR1$=ZEVGeiX=K^cCtHpAo7DI~|3n*GTJ@Q&T99QC2&1Sm@^WI^XHZ8AT49r@!MPu- z@dRKVHrxM)qxJjzNhv8;SF z%Z1_33Q7n1yj?Sc*nG7DRTs2h;Q*Y5H8I5SwW=HBq!VR`iXi6*k`8PVQG1{}7n+yt z7g0y4&`+X#nC%kkq(iI(xsVemlLoRW(djZeO7hjPIS~1B5ViUU8v6^#{#Qhevj}MM zOjL*P9ilUc;tL2DVUq)0erQFfMMLYnf5>HW(G2$R5s2G1BRBQ^5ZdDlt^ci5ItS5n zDhmEa+Zs)M2z&hYjrNDNDd|DGgsK>Rf_C46KsVAAbzKY+DO$M3fO2R7ri(Nc8cmY9 zibxVxJsz%=j-i*Ar_723j#%V@d|88oXYnQh21t_nD8uKO^bDO|dQc3%4q@LGAyNii zruoBAbOQ!NDK7cep1zXsZxepb58z>~uik^$I;fl!o3{?>PnMuf`s-pDtw?+!_c0zx zZ(IaSYD1vm1lR&7*iv8#5`M}|0_mre4fIpoxWEV~BJ!E=31~U7JRK62aWliQ1u{lH zpx>EvNtxDNy@x}rTWdfC*g>)f z`+9Ovn9}v{1{bjgL4rTHqC|BkuoPq-rNyx_#z}7oS7iCb7%fOo)nFPp2eK*13V@cK zmQ{6>xmq9Oy@pcv67mcLGSHS0_6U6-sj(938^NS!fH`f$9L!f>tD)O>Lw2kq90%+r6znxph!gFIKw;(_pbyXpKZHC$se(Ke z3#g1QKi)|b{RyBF`2{Lp386AHd|abMFla_|yLAC3ynkUd z0ZD4JHPFxX?Dz9buhXf~k!8^b^j}88JJtg|cOC3DQS4T)Ct}4m(3dvab3&l+oz@=k zC+K&;|7Y;f1^$n|5%?3?(hk)tbu9zfRLC-T%lc)&MSxT;j1Vcy!VvJ51(6FE2K)~2 zhbS2f1F^C&-q%TkTwl;NQm`=c-gkrxV?F8$Wjb9*n+v0Rs&9j*J5_ZgIO}jVBDyOg zs_jS7-f%|owj|Eftc2^}md^a%$RnYQcx~eVL8k-TXYAcvE;NL)_rlo#*c^h8P<8@r zsoP^-LL6)#rWDlm0o+X_&Jy5QRX3JZAUCug=yu-$=-LC<18$^tsDct|rx)BP{sSXP z3<{$fMGG`wp|AkEH}42QN2nl0?|504r;>nsNdbutFZooYF#zyLE#DV{2XfGU4nn2m zGm%F0GVB^mMLR*$k9w1<=^K%&U?~AwNdc|QGC*rZ2ZNSFUa+Pr6_^+7>IVS4u+*2( zWz<$mc?TYcbpxmaZd|pwTD4Cn$dfbbIGivsUW(bt(c6E5A@C-sBnBNIB$~5SrHupa zWtJzs_`9rS=eRel=y}>uuo#t%$EQo#E?o0@RHXEOGX1~7N6`sV9T*Jq4b#v@lYf!A z0{tRtR`m8d!cS~!2(mZ;#W#~q41kpSX$yyvh>8`kI%?J+a(SJe&E(Q8(g?D0Y;Ej%6STM z9WQ@|ntDSXjzb*VO(Ka-R;Ea1hc&|9nfWC<2=YGCTMG0IK2c3*))NnxJfg&#Wn2Y%l;3JpSHru?-?m5^^HeABEh0ZP{TlV%Z_T zDIr;Q>1l6dy16z{;R>4cil;ACYG6&VsA#C=99Nr!yAkoBWD#Fh65_?w z98&{RP$*#HUJ`0!2<0d=12#=#l9J2&6gE{D6zM6RcYqx&4J365X{#{CT!4A#4NhRr z6eCDXdEPC!18xxq%o49a{O1m?4pDH2V>iQo(jB=P!6BSj6xtx~B*bkX0#+elTw-U3}Q}Tqg~JuoKW@L9nE=>Og;C3Hl>KwzuQ7SS~{c)-;3?hpkFQ(a=@9 z;#}~OA>m%yt`66UDmVdk3Mw$`v38bYTD7Fw5t`it}ajpp;O{&}AD(+&lT%Lq!8^Fe%tL>QeDgykq*`)$et4oymVGZCdfYyt_#|=%PQ_g zj23kBGOL%+RP4@z(`g}U-%AJN2uQBV#URD5v6(zc2q_@W*x61{5WW6VtL#*9W(HQZ zXdDd66Hvm!%4kCrQnZ&E<@2BC_2Sb&Y9tpJa8N1ZKg$@(~GbildZ zY+htl`_)oYW_~cmvS701VrOu0tAdT=RU2yTf;8~5D)@PhCR?71d7C%FyeeByNB42? z@e{n;0u8f9xp*)TQq(MPO#v81Tb>9p4J3OZ1A{FNXS^UKKwc&9!ilg(R&pOwxw!>~ zSU9o(jDzj^6oymIn{i{7ZI>!`PCwGYzUiE+Yna z6sG>eVT-W+;^ zDidzj9f!XeWnsW2v%SsT;)4C{%y-P(j`weMVKip$0E2BP$rDQKlz5ns^~LX!`0;h| zyCojRV}1GeAby*>K7KFaw{H-?PwIg=T3hJIeABwnRz9Ig&zf=4P zf0y{1{oVADYnypDd;|VL_&ebrg1-y?B>dg*52K!)tGoQIH=-ZkzOqNzZTs6+ci)IS z!Zz>q_pI(QkNSH_J^2X1aB{D|P0?T`2xUD&Dyfk^5_ws34qe9V^ZN{b7?)#&2B9_X z4QY+N9jO-7F4&KiO05DjC|WVVOjYhAPzOH$D=^lr{Ch<-HZud_05SyilI-_d6mv*=#X>p*ZJJ(0mfQ=iG zBKA+w*+gd)j<>aa)QztQ#h#St;?A1MNpBFw;?lb{VruBDIy(R*h-waq|IyQpCS-&V z&)YB#=}2Rs4?H$!gWQdH<7Ys!Y?x;RCrZ|r+;uvv%WdK$Yj?xAd7%j2y6Z7@j2tm{ z$O9z+j+);o&m?-(!etaxl2*aji@8cqSD=HKo5Mjo+!ccNPPX5Ry|uqwID&oo!#2Fh z{3D@^u?<=PfsE98)7A_6P= zZzR6qJ+RR>-`}x^glcn6Os>_Wh~+(yI@Uk}ng&Ty_PVkZUdXh@rk3_~iERv0A;M99 zfK3%<@-X6pB#e_PE;N{2hCvo~gRsrSVQL%>{I|?`fsQO5q4Jbu;=Q-@GMUU>@-qv4A*e2290bzePJ&=SHyBwoikca&4#6mO(a-%cnO z!?MzlFJ9EP|9vj&-APQ`?s*S6J%$lT<97x>=VQr^7(<(nRA3>969wv+G!(v}$QZ;F z$t@7Mv^rAplv3@c1(G7Ek;V%cA`3Xi1mgJt5YI^^!f;bUJi~%CB%affi^TJ3h88kH zJSRY|;Yn|e4b9t!^i1ST+x*i z8`Yv%tqxv?4U|Z~x?_hBn_d??jM&V&*b&5L*Trr^Y;IlbC}Q*LV#g3$U~E^r#;@xj z(6xoUKR~kf8pS&)bBjL!wqUEYe4Anp=0U)7WL66#YM$E-f*i!NTg3Q%373lD{47L2 z+izgbqLCMV&rRbDEywn6Roa4BX}$1!d;&o|c(t8x#ND(C@>{Ocz~y8wNA1tx{W~zb zIE8^Wge0!a&=^kc?gS~n3xxQPUqmZD4$^+eC(|+FrQPv`34g@hZ0^2>%gVr74S}`V z@^WtWz6H6r zgR#(%yMo;NdULxw77X4$DE4eS7Iyk0yoFH4bO*>ooDjk2?pWPz-W`Va$V{klKFGSv zL)Q~w+Flc6vgF@)gM{Gz)dM&fSVS#6XS`#fBy)7w--J^4T-VJb=F#v?cQHR(-_X%U zAs^SZERWj`J1sNsyJO4lzpme=P3=}Za9!(q|J~Bw{;^15EZ%~f7mr_q@&)bL31%!u zCT$`jDG#<&R(F{X-LX9nzhQe$bhl^Wpn1|frRHQ++U7mDZamW6Gi9E>s^8|M4%z?twVS*x^8ZB{wRb61-dz1z|Ry3{B1Cl1n+oKvvR*#Y$w2fwcb<-<-hdlc3(vm)o9vHcrJjL+K{jh|SBlZ5IjT1l4hb5_dfBvFsj1VdW4AJd z=dhETN`?|egf)d@fiXSJ?lGryERWKMjHxqWvAa9XpeUd+FXsh}M>f(hz z&&OUiSh2=49C#5&{W`OJ7n~qP11?-D?1*|k$kZUaf&vAijv#)1*#(dFF=m63rwaKE zt{QQ1WmugM%iF79i;erD!~#ONY4XreR5|Jz5;16ekx_cXtSt&(4S0s-w$iM@a34ur zHueCUR-`jXXNV5&$Bx(@i^rU$^Y%C%gLoATp=Cqfk2ck)Zq5V*`%~n znX}|NFy7GSTs+b*LFt?|6$JawC`7pj4W(u&19Ok9Vk1r76gQZ*IZ)w4N1@GL`{FKB zTz;a%&A@a4x+z$)6<|a{Lk5=9;PZAh^1Znj?6P50LM@ceyHJro!#AN`>ZNceF|TvlYYq%DRp@zxk9AlpWu&zLD4+tixqs zdl+F#Jv)r>|J1wnN*#Nth&^CJYc!9WU7$VU@kkz??i@>XJ@W===!cLF+iDyxf7GX?wR+~U`}a_j1kDqU zUS8F3geM6y6^tJAF*tN1Htc!PDM5z>#w1U&LLJaV9b}zxL43$$4*JP9GbRUta2Q93 z)(q}E6wMLugoZ5x^+&il8v2(6h(M-hsaNWxLN~LPx|$E!rp#)=gdT$>RtHQRj+apJ z1_{Z`DlT{mV{kY~YK6T5+Pf9R6E2#U7G_pQ!Zf)7rmJl`@IK?VZR<|$Z`;oBt=(}KadSD~QN2ePX0SLu z<8dRV#3O$@ULD}-nxoQ&7M%e*$Px!W-YZ&j-YeLLRi+k;_f2lG=m zg|RjeeoGkM3|k$DJXUv_`@@{OI_pQV=yJHOYi=HvzTaz}xJITmTtR8FZjXr;>kMr2 zroB&`fruWLZ5i(6E+a3TsunF3tVnP|5;r`GPC+)d$&LX)c-Ys(4JpDOP&Gn(izDJ5 zHJ)XmBhM-S=$akOBkHaK+)yJK9^Oyy=zevn)lv0$WdGWGz`|r4|Q^|%C{q5F*$e)}H66*FCFQbNm+{UfwaTo%FVxO^> zvEO=weS8iWoUgGyau4bfvlmG}zzj+q>xQ8a$#453aDp^VD&)T$?ETD7lGT0*Dc;~T z|2l&Au$o55b&M*CRBqKk9OnfG4I3R2w+X!y4s>-G4*M$K^C?t$^PjWc?_j-S-9&PC zErW{g?y3aL5KS9}+O|K)!g)jZ3zD>`#{MUWgSnCx9{uvR_R!k2i3yu1O%jIlbEq)< zAB^9@sXg1}3kK*&*FVY@U*oHWH=@PIx@hroS4-udI9w*dy;(Ikju0YVX9QMD#QrY~ z4Z{ilp5l*A!roK(AY4O-!|H)vMz?rKiRVcGwNDL)2hx(ZsM^C9T5rJ_=M1j%BZLOYE$JC!KnG~8^g ze28FALG4D$0lrWW5LLiN#3y}$81b1G7TTyME+LeHonHhOKH_{9Yzn#S;LtB(a2X@g zV6%rBlEe7V3Qz~IT|qq=Tmc&!7Ml?j0~dCPH+W5~iqFCL6qh(|!L6sbLK1h6;fH#p zr5XAr>5>BKE`)W%Mg){V0XqhS|KUB$I5aHB4sV9Fy_l#$Cph_j?J{gg6|M zB-0|QJ^VwAgL4N^$&uHc50bdIbs5Ej7_cOnkVX4RHs~%Ec$%T7=!lhIfuTJ}C>7=W zWDs3e>*2FZ_+=!3ZEDWrLR`4ZLG-)?DnS(23NOPP*}~w48bgxe~T~< z24FKt`1ao%Xdu}WUoED_NxEy^LmTOEJJ=O_!DL5)yS^bo3|+z^6P&;?z`9(nG`yH` zx8dQB2$X2#_Wbr5OKksBl=gNU2D*d<$1;_~s?E={_rlYUJp#NhqN*NFd9t*-pEVX{ z@BwMk0kufh3&9^0C8Xf zc~b_@^}?tT!SKkhgw?TM8)mXVJdkg>7*loaXBHOpv^8nqXVU~Z(67{Mi-EpOSL#$p zhV#I6v|XrPbjtXXozX6y_s>}@F59*8js_XojulHP;iD?df8Kx$0RLtJGmRt41ME?0 zizj+VgP0sH=D!&Q8uDltj@m!Oh6`(ZG{|GC2&V^Pk)W{lNNEWsA1t~3s~n|1Yy1J2 z33db33E)}hATnm+P@f6|I}X+lOe?Md&uL&<5hipSU|ykY)pgkID`hRzwJ^oh)4k7` zG+A6gK7J#pF&F#Cmf)X+@On3@gU1tOrTd?u-bLJ6?WzSR^Hn{iN=nd$EWER!oMBib#kJDhH)wXwghZ6*EH5H1O&p?EA=$g+D~n8?Oflm(yX{FvyB@U zGv)GP(_F%KSPpXK@(Wlwr2KF@;_m=x<({i z6sY+htmQ17{|txL_|9GD5dKeV7x8(vlo7|DfNk`P5|!>R8hlIA)q$Li|Y1kiR^F?OXU72i-j*XQC`grXn zFivUJVT@8eEB9scj#=K{6KMEDq!|9}D@ar7g@631v}Di2-yCt5b#OfczygNS5ZnoJ z%E=^GyG(I@ zl8++B{4Ion*cl8Bfx`YWAt;HFg7tE%3)u(&5q5&qLkjpgh%knLgXI)ghwtSWG|^7y zOX(H1&lr|s+g-a}IF|ZyYReE13Ky#WI%|Q&H;@$UzgYqZLx=rH>b +#include "TApplication.h" // ROOT app loop +#include "ClassTransfer.h" // Reaction kinematics and MC event generation +#include "ClassAnasen.h" // ANASEN detector model classes (SX3, PW, etc.) +#include +#include +#include +#include "TLegend.h" +#include "TH1D.h" +#include "TObjArray.h" +#include "TBranch.h" +#include +#include + +//======== Generate light particle based on reaction +// calculate real and reconstructed tracks and Q-value uncertainty + +// Function to load energy loss table from file +TGraph* LoadELoss(const char* filename) { + TGraph* g = new TGraph(filename, "%lg %lg"); + return g; +} + +bool IsDeadAnode(int id){ + static std::set dead = {}; // add dead anode IDs here, 0-23 + return dead.count(id); +} + +bool IsDeadCathode(int id){ + static std::set dead = {}; // add dead cathode IDs here, 0-23 + return dead.count(id); +} + +bool IsDeadSX3(int id){ + static std::set dead = {}; // add dead SX3 IDs here, 0-23 1,7,9,3 + return dead.count(id); +} + +static std::set> ReactionProductb = { {1,1} }; // add reaction product b (light particle) A,Z pairs here, e.g. {1,1} for proton, {4,2} for alpha + +int main(int argc, char **argv){ + + printf("=========================================\n"); + printf("=== ANASEN Monte Carlo ===\n"); + printf("=========================================\n"); + + // number of events can be overridden from command line + int numEvent = 1000000; + if( argc >= 2 ) numEvent = atoi(argv[1]); + + // load energy loss tables (assume units: E in MeV, dE/dx in MeV/(mg/cm^2), density in mg/cm^3) + TGraph* elossAlpha = LoadELoss("../ELoss/E_vs_x_alpha.dat"); // for light particle (alpha) + TGraph* elossProton = LoadELoss("../ELoss/E_vs_x_proton.dat"); // for heavy particle (proton) + TGraph *invgAlpha = new TGraph(elossAlpha->GetN(), elossAlpha->GetY(), elossAlpha->GetX()); + TGraph *invgProton = new TGraph(elossProton->GetN(), elossProton->GetY(), elossProton->GetX()); + + + //Plot energy loss tables (sanity check), vis will not work if this is ran without X11 display (e.g. on cluster), so comment out if running in headless mode + auto c1 = new TCanvas("c1", "Graph Example", 800, 600); + auto g = elossAlpha; + g->SetTitle("Energy Loss Table (Alpha);cm;Kinetic Energy (MeV)"); + g->Draw("ALP"); + g->SetLineColor(kRed); + //c1->SetLogy(); + //c1->SetLogx(); + c1->Print("eloss_alpha.png"); + + auto c2 = new TCanvas("c2", "Graph Example", 800, 600); + auto g2 = elossProton; + g2->SetTitle("Energy Loss Table (Proton);cm;Kinetic Energy (MeV)"); + g2->Draw("ALP"); + g2->SetLineColor(kBlue); + c2->Print("eloss_proton.png"); + + // Reaction setup: projectile + target configuration, energy, and product IDs + TransferReaction transfer; + + transfer.SetA(18, 10, 0); // e.g., 24Mg (Z=12) with 0 excitation + transfer.SetIncidentEnergyAngle((4.4), 0, 0); // arguments are KEA in MeV/u, theta and phi in degree + transfer.Seta( 4, 2); // identify reaction product a in internal indexing e.g., 4He (alpha) + transfer.Setb(ReactionProductb.begin()->first, ReactionProductb.begin()->second); // identify reaction product b e.g., 1H (proton) + transfer.SetB(21, 11); // identify reaction product B e.g., 23Na (Z=11) + + // TODO add alpha source or alternative reaction channel selection + + // Excited state lists (target and projectile/excited products) + std::vector ExAList = {0}; // projectile excitation states in MeV + std::vector ExList = {0}; // target excitation states in MeV + + // define vertex position uniform distribution ranges (mm) + double vertexXRange[2] = { -5, 5}; // mm + double vertexYRange[2] = { -5, 5}; + double vertexZRange[2] = { -100, 100}; + + // detector resolution / uncertainty parameters + double sigmaSX3_W = -1; // mm, if < 0 use mid-point (no spread in SX3 horizontal dimension) + double sigmaSX3_L = 3; // mm, vertical spread for SX3 + double sigmaPW_A = 0; // normalized anode uncertainty term (0-1) + double sigmaPW_C = 0; // normalized cathode uncertainty term (0-1) + + // status printout + printf("------------ Vertex :\n"); + printf("X : %7.2f - %7.2f mm\n", vertexXRange[0], vertexXRange[1]); + printf("Y : %7.2f - %7.2f mm\n", vertexYRange[0], vertexYRange[1]); + printf("Z : %7.2f - %7.2f mm\n", vertexZRange[0], vertexZRange[1]); + printf("------------ Uncertainty :\n"); + printf(" SX3 horizontal : %.1f\n", sigmaSX3_W); + printf(" SX3 vertical : %.1f\n", sigmaSX3_L); + printf(" Anode : %.1f mm\n", sigmaPW_A); + printf(" Cathode : %.1f mm\n", sigmaPW_C); + printf(" num_eve : %d \n",numEvent); + + // calculates energy/momentum/kinematics constants for transfer reaction + transfer.CalReactionConstant(); + + int nExA = ExAList.size(); + int nEx = ExList.size(); + + // optional visualization control: pass "vis" as 3rd arg + bool enableVis = (argc >= 3 && strcmp(argv[2], "vis") == 0); + TApplication *app = nullptr; + if(enableVis){ + app = new TApplication("anasenVis", &argc, argv); + } + + // storage for tracks during simulation (for visualization) + std::vector visTrackVertex, visTrackDir, visTrackHitPos; + std::vector> visTrackWires; // {anodeID, cathodeID} + + // create detector representation in memory + ANASEN * anasen = new ANASEN(); // top-level detector object + SX3 * sx3 = anasen->GetSX3(); // silicon array part + PW * pw = anasen->GetPW(); // proportional wire chamber part + + // output file + tree + TString saveFileName = "SimAnasen1.root"; + printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data()); + TFile * saveFile = new TFile(saveFileName, "recreate"); + TTree * tree = new TTree("tree", "tree"); + + + // beam and CM variables saved in tree + double KEA; + tree->Branch("beamKEA", &KEA, "beamKEA/D"); + + double thetaCM, phiCM; + tree->Branch("thetaCM", &thetaCM, "thetaCM/D"); + tree->Branch("phiCM", &phiCM, "phiCM/D"); + + // outgoing particles in lab frame (light/heavy) + double thetab, phib, Tb; + double thetaB, phiB, TB; + double dEb; + double dEB; + std::array T; + tree->Branch("thetab", &thetab, "thetab/D"); // polar angle of light particle in lab frame + tree->Branch("phib", &phib, "phib/D"); // azimuthal angle of light particle in lab frame + tree->Branch("Tb", &Tb, "Tb/D"); // kinetic energy of light particle at vertex (before energy loss) + tree->Branch("thetaB", &thetaB, "thetaB/D"); + tree->Branch("phiB", &phiB, "phiB/D"); + tree->Branch("TB", &TB, "TB/D"); // kinetic energy of heavy particle at vertex (before energy loss) + tree->Branch("dEb", &dEb, "dEb/D"); + tree->Branch("dEB", &dEB, "dEB/D"); // placeholder for heavy particle energy loss, currently set equal to light particle loss for simplicity + tree->Branch("T", &T, "T/D"); // placeholder for true Q-value, currently set to 0 for simplicity + + // excitation state identifiers + int ExAID; + double ExA; + tree->Branch("ExAID", &ExAID, "ExAID/I"); // projectile excitation state ID + tree->Branch("ExA", &ExA, "ExA/D"); // projectile excitation energy in MeV + + int ExID; + double Ex; + tree->Branch("ExID", &ExID, "ExID/I"); // target excitation state ID + tree->Branch("Ex", &Ex, "Ex/D"); // target excitation energy in MeV + + // true vertex position in target volume + double vertexX, vertexY, vertexZ; + tree->Branch("vX", &vertexX, "VertexX/D"); // true vertex X position in mm + tree->Branch("vY", &vertexY, "VertexY/D"); // true vertex Y position in mm + tree->Branch("vZ", &vertexZ, "VertexZ/D"); // true vertex Z position in mm + + // reconstructed SX3 hit position + double sx3X, sx3Y, sx3Z; + tree->Branch("sx3X", &sx3X, "sx3X/D"); // reconstructed X position from SX3 (with optional smearing) + tree->Branch("sx3Y", &sx3Y, "sx3Y/D"); // reconstructed Y position from SX3 (with optional smearing) + tree->Branch("sx3Z", &sx3Z, "sx3Z/D"); // reconstructed Z position from SX3 (with optional smearing) + + // PW nearest and next nearest wires + int anodeID[2], cathodeID[2]; + tree->Branch("aID", anodeID, "anodeID/I"); // anodeID[0] is nearest anode wire, anodeID[1] is next nearest anode wire + tree->Branch("cID", cathodeID, "cathodeID/I"); // cathodeID[0] is nearest cathode wire, cathodeID[1] is next nearest cathode wire + + // distances to nearest wires + double anodeDist[2], cathodeDist[2]; + tree->Branch("aDist", anodeDist, "anodeDist/D"); + tree->Branch("cDist", cathodeDist, "cathodeDist/D"); + + // SX3 channel assignment and Z fraction (depth) information + int sx3ID, sx3Up, sx3Dn, sx3Bk; + double sx3ZFrac; + tree->Branch("sx3ID", &sx3ID, "sx3ID/I"); + tree->Branch("sx3Up", &sx3Up, "sx3Up/I"); + tree->Branch("sx3Dn", &sx3Dn, "sx3Dn/I"); + tree->Branch("sx3Bk", &sx3Bk, "sx3Bk/I"); + tree->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D"); + + // reconstructed angles from PW track fit, method 1 and 2 + double reTheta, rePhi; + tree->Branch("reTheta", &reTheta, "reconstucted_theta/D"); + tree->Branch("rePhi", &rePhi, "reconstucted_phi/D"); + + double reTheta1, rePhi1; + tree->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D"); + tree->Branch("rePhi1", &rePhi1, "reconstucted_phi1/D"); + + // reconstructed vertex Z from PW fit + double z0; + tree->Branch("z0", &z0, "reconstucted_Z/D"); + + TTree* tree2 = tree->CloneTree(0); + tree2->SetName("tree2"); + + //========timer + TBenchmark clock; + bool shown ; + clock.Reset(); + clock.Start("timer"); + shown = false; + int ELossTotal = 0; + + //================================= Calculate event loop + for( int i = 0; i < numEvent ; i++){ + + // randomly sample target/projectile excitations + ExAID = gRandom->Integer(nExA); + ExA = ExAList[ExAID]; + transfer.SetExA(ExA); + + ExID = gRandom->Integer(nEx); + Ex = ExList[ExID]; + transfer.SetExB(Ex); + + // recalc kinematic constants for chosen states + transfer.CalReactionConstant(); + + // isotropic CM direction + thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ; + phiCM = (gRandom->Rndm() - 0.5) * TMath::TwoPi(); + + //==== Calculate reaction kinematics in lab frame + TLorentzVector * output = transfer.Event(thetaCM, phiCM); // returns array of outputs + TLorentzVector Pb = output[2]; // light particle or product A + TLorentzVector PB = output[3]; // heavy particle or product B + + thetab = Pb.Theta() * TMath::RadToDeg(); + thetaB = PB.Theta() * TMath::RadToDeg(); + + Tb = (Pb.E() - Pb.M()); // kinetic energy of light particle at vertex (before energy loss) units of MeV + TB = (PB.E() - PB.M()); + T[0] = Tb; + T[1] = TB; + //if (Tb < 1.5) { + // //skip event if light particle energy after loss is below detection threshold of 1.5 MeV + // continue; + //} + + phib = Pb.Phi() * TMath::RadToDeg(); + phiB = PB.Phi() * TMath::RadToDeg(); + + // vertex position in target volume + vertexX = (vertexXRange[1]- vertexXRange[0])*gRandom->Rndm() + vertexXRange[0]; + vertexY = (vertexYRange[1]- vertexYRange[0])*gRandom->Rndm() + vertexYRange[0]; + vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() + vertexZRange[0]; + + TVector3 vertex(vertexX, vertexY, vertexZ); + + // set direction vector from lab angle + TVector3 dir(1, 0, 0); + dir.SetTheta(thetab * TMath::DegToRad()); + dir.SetPhi(phib * TMath::DegToRad()); + + // run detector response models for PW and SX3 + pw->FindWireID(vertex, dir, false); + sx3->FindSX3Pos(vertex, dir, false); + + PWHitInfo hitInfo = pw->GetHitInfo(); + + anodeID[0] = hitInfo.nearestWire.first; // nearest anode wire ID + cathodeID[0] = hitInfo.nearestWire.second; // nearest cathode wire ID + anodeID[1] = hitInfo.nextNearestWire.first; // next nearest anode wire ID + cathodeID[1] = hitInfo.nextNearestWire.second; // next nearest cathode wire ID + + anodeDist[1] = hitInfo.nextNearestDist.first; // distance to next nearest anode wire + cathodeDist[1] = hitInfo.nextNearestDist.second; // distance to next nearest cathode wire + + if(IsDeadAnode(anodeID[0])) continue; + if(IsDeadCathode(cathodeID[0])) continue; + + // SX3 hit channel info and depth fraction + sx3ID = sx3->GetID(); + + if(IsDeadSX3(sx3ID)) continue; + + anodeDist[0] = hitInfo.nearestDist.first; // distance to nearest anode wire + cathodeDist[0] = hitInfo.nearestDist.second; // distance to nearest cathode wire + + if( sx3ID >= 0 ){ + sx3Up = sx3->GetChUp(); + sx3Dn = sx3->GetChDn(); + sx3Bk = sx3->GetChBk(); + sx3ZFrac = sx3->GetZFrac(); + + // apply intrinsic detector resolution to true SX3 hit position + // for no smearing comment out and use GetHitPos(); + TVector3 hitPos = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L); + + sx3X = hitPos.X(); + sx3Y = hitPos.Y(); + sx3Z = hitPos.Z(); + + // store track data for visualization if enabled + if(enableVis){ + visTrackVertex.push_back(vertex); + visTrackDir.push_back(dir); + visTrackHitPos.push_back(hitPos); + visTrackWires.push_back({anodeID[0], cathodeID[0]}); + } + // reconstruct track from PW readings + SX3 hit + pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false); + reTheta = pw->GetTrackTheta() * TMath::RadToDeg(); + rePhi = pw->GetTrackPhi() * TMath::RadToDeg(); + + // alternative track algorithm with uncertainty parameters + pw->CalTrack2(hitPos, hitInfo, sigmaPW_A, sigmaPW_C, false); + reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg(); + rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg(); + + z0 = pw->GetZ0(); + dEb = 0; + dEB = 0; + tree->Fill(); + + //Energy loss + double dl = (hitPos - vertex).Mag(); // path length in units of cm + if (numEvent <= 100){ + //printf("Event %d: Ekin before loss = %f MeV, distance = %f cm\n", i, Tb, dl); + //printf("Total T before loss: %f MeV\n", T); + } + double tb_temp = Tb; + + dEb = tb_temp - Tb; // total energy loss + if (ReactionProductb.count({4, 2})){ // if light particle is alpha, use alpha energy loss table + double x0b = invgAlpha->Eval(Tb); + x0b = x0b + dl; + Tb = elossAlpha->Eval(x0b); + } else if (ReactionProductb.count({1, 1})){ // if light particle is proton, use proton energy loss table + double x0b = invgProton->Eval(Tb); + x0b = x0b + dl; + Tb = elossProton->Eval(x0b); + } else { + // for other particle types, can add additional energy loss tables or use a generic approximation + // for now, we will just apply a simple linear energy loss as a placeholder + double dE_dx = 5; // MeV/cm, placeholder value for energy loss per unit length + Tb = Tb - dE_dx * dl; + } + + //if (Tb < 0) { + // Tb = TMath::QuietNaN(); + //} + + dEb = tb_temp - Tb; // total energy loss + + // fill tree2 with energy loss adjusted data + //Fill T so it can make a histogram of both Tb and TB in root script + T[0] = Tb; + T[1] = 0; + //to plot both as one histogram in root, can use tree2->Draw("T(0)"); for light particle and tree2->Draw("T(1)") for heavy particle + + tree2->Fill(); + + if (numEvent <= 10){ + //printf("Event %d: Tb after energy loss = %f MeV, energy loss = %f MeV\n", i, Tb, tb_temp - Tb); + } //to give in scientific notation, use %e instead of %f in the printf format string. For example: printf("Event %d: Tb after energy loss = %e MeV, energy loss = %e MeV\n", i, Tb, tb_temp - Tb); + ELossTotal += (tb_temp - Tb); + + }else{ + // no valid SX3 hit: mark clearly invalid + sx3Up = -1; + sx3Dn = -1; + sx3Bk = -1; + sx3ZFrac = TMath::QuietNaN(); + + sx3X = TMath::QuietNaN(); + sx3Y = TMath::QuietNaN(); + sx3Z = TMath::QuietNaN(); + + reTheta = TMath::QuietNaN(); + rePhi = TMath::QuietNaN(); + reTheta1 = TMath::QuietNaN(); + rePhi1 = TMath::QuietNaN(); + z0 = TMath::QuietNaN(); + dEb = TMath::QuietNaN(); + dEB = TMath::QuietNaN(); + //Tb = -12354567; // mark kinetic energy as invalid for no hit case + // fill tree with original data (no energy loss for these events) + //comment out tree fill for no hit case + //tree->Fill(); + //tree2->Fill(); + } + + //#################################################################### Timer + // measure elapsed real time and print progress roughly every 10 sec + clock.Stop("timer"); + Double_t time = clock.GetRealTime("timer"); + clock.Start("timer"); + + if ( !shown ) { + if (fmod(time, 10) < 1 ){ + printf( "%10d[%2d%%]| %8.2f sec | expect: %5.1f min \n", i, TMath::Nint((i+1)*100./numEvent), time , numEvent*time/(i+1)/60); + shown = 1; + } + } else { + if (fmod(time, 10) > 9 ){ + shown = 0; + } + } + + } + + // write results to ROOT file and close + //tree->Write(); + //tree2->Write(); + tree->Write("", TObject::kOverwrite); + tree2->Write("", TObject::kOverwrite); + int count = tree->GetEntries(); + int count2 = tree2->GetEntries(); + saveFile->Close(); + + printf("=============== done. saved as %s. tree entries: %d, tree2 entries: %d\n", saveFileName.Data(), count, count2); + printf("Total energy loss across all events: %f MeV\n", (double)ELossTotal); + printf("Average energy loss across events: %f MeV\n", (double)ELossTotal / count); + + if(enableVis){ // to enable visualization, run with 3rd argument "vis", e.g. "./anasenMC 1000 vis" + printf("Displaying geometry with %zu tracks from simulation\n", visTrackVertex.size()); + + // Build full geometry with all wires + anasen->DrawAnasen(0, 23, 0, 23, -1, true); + + // Add all stored tracks to the geometry + TGeoManager *geom = anasen->GetGeoManager(); + TGeoVolume *worldBox = anasen->GetWorldBox(); + + if(geom && worldBox && visTrackVertex.size() > 0){ + int trackNodeID = 500; // start node IDs for tracks + + for(size_t iTrack = 0; iTrack < visTrackVertex.size(); ++iTrack){ + TVector3 vertex = visTrackVertex[iTrack]; + TVector3 dir = visTrackDir[iTrack]; + TVector3 hitPos = visTrackHitPos[iTrack]; + + double theta = dir.Theta() * TMath::RadToDeg(); + double phi = dir.Phi() * TMath::RadToDeg(); + + // Add a line marker at the vertex + TGeoVolume *startMarker = geom->MakeSphere("startMarker", 0, 0, 2.0); + startMarker->SetLineColor(kBlack); + worldBox->AddNode(startMarker, trackNodeID, + new TGeoCombiTrans(vertex.X(), vertex.Y(), vertex.Z(), + new TGeoRotation("rot", 0, 0, 0))); + trackNodeID++; + + // Add track line from vertex toward hit position + TGeoVolume *trackLine = geom->MakeTube("trackLine", 0, 0, 0.08, 150.0); + trackLine->SetLineColor(kBlue); + worldBox->AddNode(trackLine, trackNodeID, + new TGeoCombiTrans(vertex.X(), vertex.Y(), vertex.Z(), + new TGeoRotation("rotTrack", phi + 90, theta, 0))); + trackNodeID++; + + // Add hit position marker + TGeoVolume *hitMarker = geom->MakeSphere("hitMarker", 0, 0, 2.0); + hitMarker->SetLineColor(kRed); + worldBox->AddNode(hitMarker, trackNodeID, + new TGeoCombiTrans(hitPos.X(), hitPos.Y(), hitPos.Z(), + new TGeoRotation("rotHit", 0, 0, 0))); + trackNodeID++; + } + + // Redraw geometry with all tracks + geom->CloseGeometry(); + geom->SetVisLevel(4); + worldBox->Draw("ogle"); + } + + if(app){ + printf("Entering ROOT event loop\n"); + app->Run(); + } + } + + delete anasen; + delete elossAlpha; + delete elossProton; + + + return 0; + +} \ No newline at end of file