python structure fixes

This commit is contained in:
James Szalkie 2026-06-05 14:14:21 -04:00
parent 92d831678e
commit e4a490245c
5 changed files with 1194 additions and 586 deletions

View File

@ -57,8 +57,8 @@ private:
const int numDet = 12; const int numDet = 12;
const float radius = 88; const float radius = 88;
const float width = 40; const float width = 40;
const float length = 174.3; // 75 const float length = 75; // 75
const float gap = 0; // 46 const float gap = 46; // 46
short id; // -1 when no hit short id; // -1 when no hit
short chUp; short chUp;

View File

@ -97,13 +97,13 @@ int main(int argc, char **argv){
if( argc >= 2 ) numEvent = atoi(argv[1]); if( argc >= 2 ) numEvent = atoi(argv[1]);
TransferReaction transfer; TransferReaction transfer;
transfer.SetA(27, 13, 0); // 18Ne projectile transfer.SetA(1, 1, 0); // 18Ne projectile
//transfer.SetIncidentEnergyAngle(0, 0, 0); // KEA in MeV/u, theta and phi in rad //transfer.SetIncidentEnergyAngle(0, 0, 0); // KEA in MeV/u, theta and phi in rad
TGraph* elossBeam = LoadELoss("../ELoss/E_vs_x_Ne-18.dat"); TGraph* elossBeam = LoadELoss("../ELoss/E_vs_x_proton.dat");
transfer.Seta(4, 2); // 4He target transfer.Seta(4, 2); // 4He target
transfer.Setb(1, 1); // outgoing proton from the primary transfer transfer.Setb(1, 1); // outgoing proton from the primary transfer
transfer.SetB(30, 14); // 21Na* heavy product transfer.SetB(4, 2); // 21Na* heavy product
const double beamA = 27.0; // mass number of 27Al beam const double beamA = 1.0; // mass number of 27Al beam
bool enableSequentialDecay = false; // turning to false to disable sequential decay for now, can be set to true to enable bool enableSequentialDecay = false; // turning to false to disable sequential decay for now, can be set to true to enable
const int decayDaughterA = 20; const int decayDaughterA = 20;
@ -118,7 +118,7 @@ int main(int argc, char **argv){
// define vertex position uniform distribution ranges (mm) // define vertex position uniform distribution ranges (mm)
double vertexXRange[2] = { -5, 5}; // mm - 5, 5 double vertexXRange[2] = { -5, 5}; // mm - 5, 5
double vertexYRange[2] = { -5, 5}; // -5, 5 double vertexYRange[2] = { -5, 5}; // -5, 5
double vertexZRange[2] = { -100, 100}; // -174.3, 174.3 (full length of gas volume, centered at 0) double vertexZRange[2] = { -174.3, 174.3}; // -174.3, 174.3 (full length of gas volume, centered at 0)
const double beamEntranceZ = -174.3; //vertexZRange[0]; // mm, assumed beam entrance into the gas const double beamEntranceZ = -174.3; //vertexZRange[0]; // mm, assumed beam entrance into the gas
@ -547,10 +547,13 @@ int main(int argc, char **argv){
vertexX2 = vertexX; vertexX2 = vertexX;
vertexY2 = vertexY; vertexY2 = vertexY;
vertexZ2 = vertexZ; vertexZ2 = vertexZ;
qqqX = TMath::QuietNaN();
qqqY = TMath::QuietNaN();
qqqZ = TMath::QuietNaN();
tree2->Fill(); tree2->Fill();
}else if (qqqID >= 0){ }else if (false) {//(qqqID >= 0){
// handle QQQ hit case // handle QQQ hit case
sx3Up = -1; sx3Up = -1;
sx3Dn = -1; sx3Dn = -1;

501
ELoss/E_vs_x_Al-27.dat Normal file
View File

@ -0,0 +1,501 @@
Distance_cm Energy_MeV
-0.0 72.0
0.07516627746635594 71.85573146292585
0.1502674587514124 71.71146292585169
0.22530353863089683 71.56719438877755
0.30027451183674336 71.4229258517034
0.37518037305649476 71.27865731462926
0.45002111693280733 71.1343887775551
0.52479673806283 70.99012024048096
0.599507230997694 70.84585170340681
0.6741525902418835 70.70158316633265
0.748732810252672 70.55731462925851
0.8232478854395513 70.41304609218436
0.8976978101635771 70.26877755511022
0.9720825787368177 70.12450901803606
1.0464021854216836 69.98024048096191
1.1206566244303215 69.83597194388777
1.1948458899240002 69.69170340681362
1.268969976012413 69.54743486973948
1.3430288767530822 69.40316633266532
1.4170225861506371 69.25889779559118
1.490951098156198 69.11462925851703
1.5648144066666436 68.97036072144287
1.6386125055239398 68.82609218436873
1.7123453885144593 68.68182364729458
1.7860130493682191 68.53755511022044
1.8596154817582156 68.39328657314628
1.933152679299634 68.24901803607214
2.006624635549164 68.10474949899799
2.0800313440041944 67.96048096192384
2.153372798102072 67.8162124248497
2.226648991219346 67.67194388777554
2.299859916670933 67.5276753507014
2.3730055677093755 67.38340681362725
2.4460859375239763 67.2391382765531
2.5191010192400345 67.09486973947895
2.5920508059179634 66.9506012024048
2.6649352905524664 66.80633266533066
2.7377544660717033 66.6620641282565
2.8105083253363703 66.51779559118236
2.883196861138875 66.37352705410821
2.955820066202387 66.22925851703407
3.028377933179984 66.08498997995991
3.100870454653685 65.94072144288576
3.173297623133536 65.79645290581162
3.245659431056685 65.65218436873747
3.317955870786368 65.50791583166333
3.3901869346109943 65.36364729458917
3.4623526147430987 65.21937875751503
3.5344529033183933 65.07511022044088
3.6064877923947 64.93084168336672
3.6784572739509387 64.78657314629258
3.750361339886098 64.64230460921843
3.8221999820181236 64.49803607214429
3.893973192082893 64.35376753507013
3.965680961733064 64.209498997996
4.03732328253702 64.06523046092184
4.108900145977691 63.920961923847685
4.180411543451432 63.77669338677354
4.251857466266871 63.63242484969938
4.323237905643703 63.488156312625236
4.394552852711519 63.34388777555109
4.465802298508576 63.19961923847694
4.536986233980569 63.055350701402794
4.608104649979376 62.91108216432865
4.679157537261796 62.7668136272545
4.750144886488255 62.622545090180346
4.82106668822149 62.4782765531062
4.8919229329252385 62.33400801603205
4.96271361096288 62.189739478957904
5.0334387125960705 62.04547094188376
5.104098227983351 61.90120240480961
5.1746921471787415 61.75693386773546
5.245220460130313 61.61266533066131
5.315683156678721 61.46839679358716
5.386080226555749 61.32412825651301
5.456411659382805 61.179859719438866
5.5266774446693985 61.03559118236472
5.5968775718116 60.89132264529057
5.66701203009048 60.74705410821642
5.737080808670504 60.60278557114227
5.807083896597937 60.45851703406812
5.877021282799196 60.314248496993976
5.946892956079183 60.16997995991983
6.016698905119601 60.02571142284568
6.0864391184772355 59.881442885771534
6.156113584582216 59.73717434869738
6.2257222917362345 59.59290581162323
6.29526522811077 59.448637274549085
6.364742381745252 59.30436873747494
6.434153740545209 59.16010020040079
6.503499292280392 59.015831663326644
6.57277902458286 58.8715631262525
6.641992924945049 58.72729458917834
6.711140980717785 58.583026052104195
6.7802231791083045 58.43875751503005
6.849239507178213 58.2944889779559
6.918189951841421 58.15022044088175
6.987074499862048 58.005951903807606
7.0558931378523 57.86168336673346
7.124645852270304 57.717414829659305
7.193332629417903 57.57314629258516
7.261953455438445 57.42887775551101
7.3305083163145035 57.28460921843686
7.3989971978655795 57.140340681362716
7.467420085745765 56.99607214428857
7.5357769654413715 56.851803607214414
7.604067822268504 56.70753507014027
7.672292641370628 56.56326653306612
7.740451407716073 56.41899799599197
7.808544106095498 56.274729458917825
7.87657072111933 56.13046092184368
7.9445312372151475 55.98619238476953
8.012425638625036 55.84192384769538
8.080253909402876 55.69765531062123
8.148016033411626 55.55338677354708
8.215711994320529 55.409118236472935
8.283341775602285 55.26484969939879
8.35090536053018 55.12058116232464
8.418402732175172 54.97631262525049
8.48583387340292 54.83204408817634
8.553198766870759 54.68777555110219
8.620497395024662 54.543507014028044
8.687729740096108 54.3992384769539
8.754895784098922 54.25496993987975
8.821995508826063 54.1107014028056
8.889028895846357 53.966432865731456
8.955995926501174 53.8221643286573
9.022896581901046 53.677895791583154
9.089730842922247 53.53362725450901
9.156498690203303 53.38935871743486
9.223200104141448 53.24509018036071
9.289835064889012 53.100821643286565
9.356403552349773 52.95655310621241
9.42290554617522 52.812284569138264
9.489341025760782 52.668016032064116
9.555709970241967 52.52374749498997
9.622012358490466 52.37947895791582
9.68824816911016 52.235210420841675
9.754417380433093 52.09094188376753
9.820519970515353 51.94667334669337
9.886555917132887 51.802404809619226
9.952525197777272 51.65813627254508
10.018427789651378 51.51386773547093
10.084263669664981 51.369599198396784
10.150032814430297 51.22533066132264
10.21573520025744 51.08106212424849
10.281370803149805 50.936793587174336
10.346939598799363 50.79252505010019
10.4124415625819 50.64825651302604
10.47787666955215 50.503987975951894
10.543244894438859 50.35971943887775
10.608546211639755 50.2154509018036
10.67378059521645 50.07118236472945
10.738948018889236 49.9269138276553
10.804048456031795 49.78264529058115
10.86908187966584 49.638376753507
10.934048262455631 49.494108216432856
10.99894757670242 49.34983967935871
11.063779794338794 49.20557114228456
11.128544886922919 49.06130260521041
11.193242825632675 48.91703406813626
11.25787358125972 48.77276553106211
11.322437124203418 48.628496993987966
11.386933424464676 48.48422845691382
11.451362451639676 48.33995991983967
11.515724174913492 48.195691382765524
11.5800185630536 48.05142284569137
11.644245584403263 47.90715430861722
11.708405206874833 47.762885771543075
11.772497397942887 47.61861723446893
11.836522124637286 47.47434869739478
11.900479353536083 47.330080160320634
11.964369050758325 47.18581162324649
12.028191181956718 47.04154308617233
12.091945712310155 46.897274549098185
12.15563260651614 46.75300601202404
12.219251828783042 46.60873747494989
12.282803342822238 46.46446893787574
12.346287111840102 46.320200400801596
12.409703098529857 46.17593186372745
12.473051265063283 46.031663326653295
12.536331573082263 45.88739478957915
12.599543983690205 45.743126252505
12.662688457443284 45.59885771543085
12.725764954341537 45.454589178356706
12.788773433819797 45.31032064128256
12.851713854738469 45.166052104208404
12.914586175374119 45.02178356713426
12.977390353409929 44.87751503006011
13.040126345925941 44.73324649298596
13.10279410938915 44.588977955911815
13.1653935996434 44.44470941883767
13.22792477189911 44.30044088176352
13.290387580722806 44.15617234468937
13.35278198002645 44.01190380761522
13.4151079230566 43.86763527054107
13.477365362383358 43.723366733466925
13.539554249889097 43.57909819639278
13.601674536757013 43.43482965931863
13.663726173459446 43.29056112224448
13.72570910974599 43.14629258517033
13.787623294631373 43.00202404809618
13.849468676383143 42.857755511022035
13.91124520250909 42.71348697394789
13.972952819744453 42.56921843687374
14.034591474038892 42.42494989979959
14.096161110543205 42.28068136272544
14.157661673595797 42.13641282565129
14.219093106708915 41.992144288577144
14.280455352554606 41.847875751503
14.341748352950415 41.70360721442885
14.402972048844822 41.5593386773547
14.464126380302396 41.415070140280555
14.52521128648868 41.2708016032064
14.586226705654767 41.126533066132254
14.64717257512162 40.98226452905811
14.708048831264072 40.83799599198396
14.768855409494517 40.69372745490981
14.829592244246307 40.549458917835665
14.89025926895683 40.40519038076152
14.950856416050263 40.26092184368736
15.011383616919975 40.116653306613216
15.071840801910644 39.97238476953907
15.132227900299984 39.82811623246492
15.192544840280139 39.683847695390774
15.25279154893872 39.53957915831663
15.312967952239482 39.39531062124248
15.37307397500262 39.251042084168326
15.43310954088468 39.10677354709418
15.493074572358102 38.96250501002003
15.552968990690355 38.818236472945884
15.612792715922655 38.67396793587174
15.672545666848297 38.52969939879759
15.732227760990545 38.385430861723435
15.791838914580092 38.24116232464929
15.851379042532098 38.09689378757514
15.910848058422765 37.952625250500994
15.970245874465451 37.808356713426846
16.029572401486334 37.6640881763527
16.088827548899577 37.51981963927855
16.148011224682026 37.3755511022044
16.20712333534739 37.23128256513025
16.266163785919936 37.0870140280561
16.32513247990765 36.942745490981956
16.38402931927487 36.79847695390781
16.442854204414374 36.65420841683366
16.501607034118933 36.509939879759514
16.56028770555229 36.36567134268536
16.61889611421955 36.22140280561121
16.677432153937005 36.077134268537066
16.735895716801355 35.93286573146292
16.794286693158295 35.78859719438877
16.852604971570504 35.644328657314624
16.910850438784973 35.50006012024048
16.969022979699705 35.35579158316632
17.027122477329716 35.211523046092175
17.085148812772392 35.06725450901803
17.14310186517212 34.92298597194388
17.200981511684237 34.77871743486973
17.258787627438227 34.634448897795586
17.31652008550021 34.49018036072143
17.374178756834638 34.345911823647285
17.431763510265277 34.20164328657314
17.48927421243533 34.05737474949899
17.54671072776682 33.91310621242484
17.60407291841911 33.768837675350696
17.66136064424662 33.62456913827655
17.71857376275566 33.480300601202394
17.7757121290604 33.33603206412825
17.83277559583795 33.1917635270541
17.889764013282555 33.04749498997995
17.946677229058807 32.903226452905805
18.003515088253952 32.75895791583166
18.060277433329215 32.61468937875751
18.116964104070153 32.47042084168336
18.17357493753596 32.32615230460921
18.230109768007797 32.18188376753506
18.286568426936025 32.037615230460915
18.342950742886398 31.89334669338677
18.399256541485144 31.74907815631262
18.455485645362934 31.604809619238473
18.511637874097723 31.460541082164326
18.567720934190692 31.316272545090175
18.623743934764097 31.172004008016028
18.679708837636362 31.02773547094188
18.7356163243514 30.883466933867734
18.791466457919757 30.739198396793583
18.847259301916758 30.594929859719436
18.902994920491153 30.45066132264529
18.958673378374048 30.306392785571138
19.014294740888168 30.16212424849699
19.06985907395749 30.017855711422843
19.125366444117223 29.873587174348692
19.180816918524187 29.729318637274545
19.23621056496759 29.585050100200398
19.29154745188025 29.44078156312625
19.34682764835025 29.2965130260521
19.402051224133082 29.152244488977953
19.45721824966429 29.007975951903806
19.51232879607264 28.863707414829655
19.56738293519384 28.719438877755508
19.622380739584862 28.57517034068136
19.677322282538874 28.43090180360721
19.732207638100814 28.286633266533062
19.78703688108367 28.142364729458915
19.841810087085467 27.998096192384768
19.896527332507027 27.853827655310617
19.951188694570526 27.70955911823647
20.005794251338905 27.565290581162323
20.060344081736183 27.421022044088172
20.114838265568686 27.276753507014025
20.169276883547308 27.132484969939878
20.223660017310788 26.98821643286573
20.277987749450123 26.84394789579158
20.332260163534116 26.699679358717432
20.38647734413619 26.555410821643285
20.440639376862485 26.411142284569134
20.494746348381316 26.266873747494987
20.54879834645411 26.12260521042084
20.602795459967833 25.97833667334669
20.656737778969063 25.834068136272542
20.71062539469974 25.689799599198395
20.76445839963473 25.545531062124248
20.818236887521277 25.401262525050097
20.87196095342046 25.25699398797595
20.925630693750765 25.112725450901802
20.979246206333894 24.96845691382765
21.032807590442918 24.824188376753504
21.086314946852937 24.679919839679357
21.139768377894377 24.535651302605206
21.193167987509028 24.39138276553106
21.246513881309063 24.247114228456912
21.299806166639137 24.102845691382765
21.353044952641763 23.958577154308614
21.406230350326144 23.814308617234467
21.45936247264069 23.67004008016032
21.512441434549384 23.52577154308617
21.565467353112222 23.38150300601202
21.618440347570004 23.237234468937874
21.671360539433685 23.092965931863727
21.724228052578532 22.948697394789576
21.777043013343416 22.80442885771543
21.8298055506355 22.660160320641282
21.882515796040614 22.51589178356713
21.935173883939697 22.371623246492984
21.987779951631616 22.227354709418837
22.04033413946274 22.083086172344686
22.092836590963664 21.93881763527054
22.145287452993507 21.79454909819639
22.197686875892206 21.650280561122244
22.250035013641273 21.506012024048093
22.302332024033543 21.361743486973946
22.354578068852387 21.2174749498998
22.40677331406098 21.073206412825648
22.458917930002215 20.9289378757515
22.51101209160986 20.784669338677354
22.563055978631667 20.640400801603203
22.61504977586508 20.496132264529056
22.666993673406367 20.35186372745491
22.718887866913878 20.20759519038076
22.770732557886358 20.06332665330661
22.822527953957145 19.919058116232463
22.874274269205262 19.774789579158316
22.92597172448435 19.630521042084165
22.97762054777056 19.486252505010018
23.029220974530553 19.34198396793587
23.08077324811076 19.19771543086172
23.13227762014925 19.053446893787573
23.183734351011562 18.909178356713426
23.23514371025193 18.76490981963928
23.286505977101466 18.620641282565128
23.337832649142197 18.47637274549098
23.389143342951172 18.332104208416833
23.44044820267803 18.187835671342683
23.491749757183765 18.043567134268535
23.5430491586445 17.899298597194388
23.594347434511615 17.75503006012024
23.64564564229102 17.61076152304609
23.696944870858683 17.466492985971943
23.748246241847326 17.322224448897796
23.7995509111086 17.177955911823645
23.850860070255433 17.033687374749498
23.90217494828952 16.88941883767535
23.95349681331927 16.7451503006012
24.004826974373863 16.600881763527052
24.056166783319547 16.456613226452905
24.107517636884623 16.312344689378758
24.158880978800095 16.168076152304607
24.21025830206343 16.02380761523046
24.261651151333425 15.87953907815631
24.313061125464696 15.735270541082162
24.364489880190984 15.591002004008013
24.41593913096712 15.446733466933866
24.467410655980185 15.302464929859717
24.518906299341218 15.158196392785568
24.570427974469663 15.01392785571142
24.621977667683645 14.869659318637272
24.67355744201016 14.725390781563124
24.72516944123037 14.581122244488975
24.7768158941763 14.436853707414826
24.828499119296573 14.29258517034068
24.880221529510163 14.14831663326653
24.93198563736862 14.004048096192383
24.98379406054891 13.859779559118234
25.035649527700766 13.715511022044085
25.087554884674383 13.571242484969938
25.139513101156414 13.426973947895789
25.19152727774457 13.282705410821642
25.243600653493615 13.138436873747493
25.295736613968387 12.994168336673344
25.347938699842445 12.849899799599196
25.400210616084376 12.705631262525047
25.452556241777383 12.5613627254509
25.50497964062188 12.417094188376751
25.557485072175194 12.272825651302602
25.610077003887486 12.128557114228455
25.662760123998297 11.984288577154306
25.715539355364164 11.840020040080159
25.768419870294395 11.69575150300601
25.82140710647933 11.551482965931863
25.874506784103673 11.407214428857714
25.927724924246363 11.262945891783565
25.981067868678664 11.118677354709417
26.034542301183215 10.974408817635268
26.08815527052934 10.830140280561121
26.14191421525389 10.685871743486972
26.19582699041243 10.541603206412823
26.24990189648321 10.397334669338676
26.304147710625887 10.253066132264527
26.35857372051921 10.10879759519038
26.41318976102666 9.96452905811623
26.468006253967225 9.820260521042082
26.523034251300246 9.675991983967934
26.57828548206925 9.531723446893785
26.633742337775963 9.387454909819638
26.689338083235043 9.24318637274549
26.745028925059557 9.09891783567134
26.800816117138876 8.954649298597193
26.85670585548451 8.810380761523044
26.912706900318735 8.666112224448897
26.96882861227707 8.521843687374748
27.02508099838088 8.377575150300599
27.081474762033395 8.233306613226452
27.13802135745477 8.089038076152303
27.194733049022997 7.9447695390781545
27.25162297604574 7.800501002004006
27.30870522355517 7.656232464929858
27.365994899794575 7.51196392785571
27.42350822115415 7.367695390781562
27.481262605415292 7.223426853707413
27.53927677428101 7.079158316633265
27.597570866307027 6.934889779559117
27.65616656150753 6.790621242484969
27.715087219095473 6.646352705410821
27.774358030034815 6.502084168336672
27.83400618633711 6.357815631262524
27.894061069335237 6.2135470941883755
27.95455445952136 6.069278557114227
28.01552077095618 5.925010020040079
28.07699731375548 5.780741482965931
28.139024588755394 5.636472945891782
28.201646619170532 5.492204408817634
28.264911324916156 5.347935871743486
28.32887094629981 5.203667334669338
28.39358252504175 5.05939879759519
28.459108452110375 4.915130260521041
28.525517093726656 4.770861723446893
28.592883509188262 4.6265931863727445
28.661290277002124 4.482324649298596
28.73082844934073 4.338056112224448
28.80159865924444 4.193787575150299
28.873712410532516 4.049519038076151
28.94729358739528 3.905250501002003
29.022480229567456 3.760981963927855
29.099426630434042 3.6167134268537064
29.178305830217344 3.4724448897795583
29.259312595666188 3.3281763527054102
29.342667002978146 3.1839078156312617
29.42861877421806 3.0396392785571136
29.517452562336686 2.8953707414829655
29.609494440447484 2.751102204408817
29.705119933643836 2.606833667334669
29.804764045611353 2.4625651302605203
29.908933891300695 2.318296593186372
30.018224771400167 2.174028056112224
30.13334084505466 2.0297595190380755
30.255122020930816 1.8854909819639276
30.384579364471794 1.7412224448897793
30.522942318747475 1.596953907815631
30.67172251674337 1.4526853707414829
30.832801142369917 1.3084168336673345
31.008549911062335 1.1641482965931862
31.20199978106649 1.0198797595190379
31.41706960071087 0.8756112224448896
31.658882845241283 0.7313426853707414
31.929915336524896 0.5870741482965931
32.22943849820706 0.44280561122244483
32.559672476181774 0.29853707414829656
32.922180535599225 0.15426853707414828
33.267604963464514 0.01

File diff suppressed because it is too large Load Diff

View File

@ -20,6 +20,8 @@ import os
import periodictable as pt import periodictable as pt
import re import re
from matplotlib.colors import LinearSegmentedColormap from matplotlib.colors import LinearSegmentedColormap
from mpl_toolkits.mplot3d import Axes3D
import nbformat as nbf
# ROOT-like styling # ROOT-like styling
plt.rcParams.update({ plt.rcParams.update({
@ -291,50 +293,65 @@ def calculate_distance_tree2(vz, theta, z_max=34.86):
return np.abs((z_max - vz) / cos_theta) return np.abs((z_max - vz) / cos_theta)
def load_tree_arrays(tree, treename, max_events=None): def load_tree_arrays(tree, treename, max_events=None):
branches = ["Tb", "thetab", "sx3Z", "vX", "vY", "vZ", "sx3ID", "aID", "cID", "aDist", "cDist"] branches = tree.keys(recursive=False)
if treename == 'tree1':
branches.extend(["sx3X", "sx3Y"])
return tree.arrays(branches, library="np", entry_stop=max_events) return tree.arrays(branches, library="np", entry_stop=max_events)
def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86): def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
data = load_tree_arrays(tree, treename, max_events) data = load_tree_arrays(tree, treename, max_events)
#mask = data["thetab"] >= 0 #mask = (~np.isnan(data["sx3X"])) & (~np.isnan(data["Tb"]))
#data = {key: value[mask] for key, value in data.items()} #data = {key: value[mask] for key, value in data.items()}
Ei = data["Tb"] Ei = data["Tb"]
Einan = np.isnan(Ei).sum()
print(f"sx3 NaN entries: {Einan}")
theta = np.radians(data["thetab"]) theta = np.radians(data["thetab"])
vX = data["vX"] vX = data["vX"]
vY = data["vY"] vY = data["vY"]
vZ = data["vZ"] vZ = data["vZ"]
if treename == 'tree1': #if treename == 'tree1':
sx3X = data["sx3X"] sx3X = data["sx3X"]
sx3Y = data["sx3Y"] sx3Y = data["sx3Y"]
sx3Z = data["sx3Z"] sx3Z = data["sx3Z"]
mask = ~np.isnan(sx3X) & ~np.isnan(sx3Y) & ~np.isnan(sx3Z) mask = ~np.isnan(sx3X) & ~np.isnan(sx3Y) & ~np.isnan(sx3Z)
else: qqqX = data["qqqX"]
sx3Z = np.full_like(vZ, z_max) # Assume sx3Z is at z_max for tree2 qqqY = data["qqqY"]
mask = ~np.isnan(Ei) & ~np.isnan(theta) & ~np.isnan(vZ) qqqZ = data["qqqZ"]
qqqnan = np.isnan(qqqZ).sum()
print(f"qqq NaN entries: {qqqnan}")
print(f"Total: {Einan + qqqnan}")
mask2 = ~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] Eisx3 = Ei[mask]
theta = theta[mask] theta = theta[mask]
vX = vX[mask] vXsx3 = vX[mask]
vY = vY[mask] vYsx3 = vY[mask]
vZ = vZ[mask] vZsx3 = vZ[mask]
Eiqqq = Ei[mask2]
vXqqq = vX[mask2]
vYqqq = vY[mask2]
vZqqq = vZ[mask2]
if treename == 'tree1': #if treename == 'tree1':
sx3X = sx3X[mask] sx3X = sx3X[mask]
sx3Y = sx3Y[mask] sx3Y = sx3Y[mask]
sx3Z = sx3Z[mask] sx3Z = sx3Z[mask]
dsx3 = calculate_distance_tree1(vX, vY, vZ, sx3X, sx3Y, sx3Z) dsx3 = calculate_distance_tree1(vXsx3, vYsx3, vZsx3, sx3X, sx3Y, sx3Z)
else: qqqX = qqqX[mask2]
dsx3 = calculate_distance_tree2(vZ, theta, z_max=z_max) qqqY = qqqY[mask2]
qqqZ = qqqZ[mask2]
dqqq = calculate_distance_tree1(vXqqq, vYqqq, vZqqq, qqqX, qqqY, 128)
#else:
# dsx3 = calculate_distance_tree2(vZsx3, theta, z_max=z_max)
sin_theta = np.sin(theta) sin_theta = np.sin(theta)
sin_theta = np.where(sin_theta != 0, sin_theta, 1e-10) sin_theta = np.where(sin_theta != 0, sin_theta, 1e-10)
radii = np.array([3.7, 4.3]) radii = np.array([3.7, 4.3])
dA = (radii[0] - np.sqrt((vX/10)**2 + (vY/10)**2))/ sin_theta dA = (radii[0] - np.sqrt((vXsx3/10)**2 + (vYsx3/10)**2))/ sin_theta
dC = (radii[1] - np.sqrt((vX/10)**2 + (vY/10)**2))/ sin_theta dC = (radii[1] - np.sqrt((vXsx3/10)**2 + (vYsx3/10)**2))/ sin_theta
# Filter out unphysical distances (negative or unreasonably small) # Filter out unphysical distances (negative or unreasonably small)
# These typically occur at large angles where trajectory doesn't properly intersect proportional counters # These typically occur at large angles where trajectory doesn't properly intersect proportional counters
@ -351,31 +368,22 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
else: else:
aID_valid = aID >= 0 aID_valid = aID >= 0
cID_valid = cID >= 0 cID_valid = cID >= 0
distance_mask = (dA > 0.1) & (dC > 0.1) & (sx3ID >= 0) & aID_valid & cID_valid
Ei = Ei[distance_mask]
theta = theta[distance_mask]
dA = dA[distance_mask]
dC = dC[distance_mask]
dsx3 = dsx3[distance_mask]
vX = vX[distance_mask]
vY = vY[distance_mask]
vZ = vZ[distance_mask]
print(f"Computing energies for {particle} ({treename})...") print(f"Computing energies for {particle} ({treename})...")
print(f" Retained {np.sum(distance_mask)} / {len(distance_mask)} events after distance filter") #print(f" Retained {np.sum(distance_mask)} / {len(distance_mask)} events after distance filter")
EA = energy_loss(particle, Ei, dA) EA = energy_loss(particle, Eisx3, dA)
EC = energy_loss(particle, Ei, dC) EC = energy_loss(particle, Eisx3, dC)
Esx3 = energy_loss(particle, Ei, dsx3) Esx3 = energy_loss(particle, Eisx3, dsx3)
Eqqq = energy_loss(particle, Eiqqq, dqqq)
Elost = Ei - Esx3 Elost = Eisx3 - Esx3
Eprop = EA - EC Eprop = EA - EC
return { return {
"particle": particle, "particle": particle,
"Ei": Ei, "Ei": Eisx3,
"thetabi": data["thetab"],
"sx3Z": sx3Z, "sx3Z": sx3Z,
"EA": EA, "EA": EA,
"EC": EC, "EC": EC,
@ -384,7 +392,13 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
"Eprop": Eprop, "Eprop": Eprop,
"dA": dA, "dA": dA,
"dC": dC, "dC": dC,
"thetab": np.degrees(theta) "thetab": np.degrees(theta),
"sx3X": sx3X,
"sx3Y": sx3Y,
"qqqX": qqqX,
"qqqY": qqqY,
"qqqZ": qqqZ,
"Eqqq": Eqqq
} }
def process_file(filename, treename, particle=None, max_events=None): def process_file(filename, treename, particle=None, max_events=None):
@ -970,6 +984,9 @@ class MyInteractiveApp(cmd.Cmd):
data = prepare_tree_data(self.tree, treename, particle, max_events=max_events) data = prepare_tree_data(self.tree, treename, particle, max_events=max_events)
Ei = data["Ei"] Ei = data["Ei"]
thetab = data["thetab"]
sx3X = data["sx3X"]
sx3Y = data["sx3Y"]
sx3Z = data["sx3Z"] sx3Z = data["sx3Z"]
EA = data["EA"] EA = data["EA"]
EC = data["EC"] EC = data["EC"]
@ -978,6 +995,10 @@ class MyInteractiveApp(cmd.Cmd):
Eprop = data["Eprop"] Eprop = data["Eprop"]
dA = data["dA"] dA = data["dA"]
dC = data["dC"] dC = data["dC"]
qqqX = data["qqqX"]
qqqY = data["qqqY"]
qqqZ = data["qqqZ"]
qqqE = data["Eqqq"]
update_plot_data(f"{particle}_{treename}_Ei", Ei) update_plot_data(f"{particle}_{treename}_Ei", Ei)
update_plot_data(f"{particle}_{treename}_sx3Z", sx3Z) update_plot_data(f"{particle}_{treename}_sx3Z", sx3Z)
@ -991,6 +1012,110 @@ class MyInteractiveApp(cmd.Cmd):
os.makedirs(base, exist_ok=True) os.makedirs(base, exist_ok=True)
print(f"Saving plots to folder: {base} ({treename})") print(f"Saving plots to folder: {base} ({treename})")
# --- clean data ---
x = np.asarray(sx3X, dtype=float)
y = np.asarray(sx3Y, dtype=float)
z = np.asarray(sx3Z, dtype=float)
c = np.asarray(Esx3, dtype=float)
qx = np.asarray(qqqX, dtype=float)
qy = np.asarray(qqqY, dtype=float)
qz = np.asarray(qqqZ, dtype=float)
qc = np.asarray(qqqE, dtype=float)
if x.size > 0 and False:
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection='3d')
sc = ax.scatter(
x, y, z,
c=c,
cmap='viridis',
s=3,
alpha=0.8
)
ax.set_xlabel("SX3 X")
ax.set_ylabel("SX3 Y")
ax.set_zlabel("SX3 Z")
ax.set_title(f"{particle} ({treename}) SX3 3D Hit Map colored by Esx3")
cb = plt.colorbar(sc, ax=ax, pad=0.1)
cb.set_label("Esx3 (MeV)")
plt.tight_layout()
out_file = f"{base}/sx3_3D_energy.png"
plt.savefig(out_file, dpi=300)
plt.show()
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection='3d')
sc = ax.scatter(
qx, qy, qz,
c=qc,
cmap='viridis',
s=3,
alpha=0.8
)
ax.set_xlabel("SX3 X")
ax.set_ylabel("SX3 Y")
ax.set_zlabel("SX3 Z")
ax.set_title(f"{particle} ({treename}) QQQ 3D Hit Map colored by qqqTb")
cb = plt.colorbar(sc, ax=ax, pad=0.1)
cb.set_label("qqqTb (MeV)")
plt.tight_layout()
out_file = f"{base}/QQQ_3D_energy.png"
plt.savefig(out_file, dpi=300)
plt.show()
else:
print("No valid SX3 points for 3D plot.")
mask1 = ~np.isnan(Ei) & ~np.isnan(thetab)
plt.figure(figsize=(7,6))
plt.hist2d(thetab[mask1], Ei[mask1], bins=200)
plt.xlabel("thetab")
plt.ylabel("Event Energy (MeV)")
plt.title(f"{particle} ({treename}) Energy vs Theta")
plt.colorbar(label="Counts")
#plt.xlim(0,30)
#plt.ylim(0,0.45)
plt.tight_layout()
plt.savefig(f"{base}/E_vs_theta.png", dpi=300)
plt.show()
mask1 = (Esx3 > 0) & ~np.isnan(thetab)
plt.figure(figsize=(7,6))
plt.hist2d(thetab[mask1], Esx3[mask1], bins=200)
plt.xlabel("thetab")
plt.ylabel("Esx3")
plt.title(f"{particle} ({treename}) sx3 Energy vs Theta")
plt.colorbar(label="Counts")
#plt.xlim(0,30)
#plt.ylim(0,0.45)
plt.tight_layout()
plt.savefig(f"{base}/sx3E_vs_theta.png", dpi=300)
plt.show()
mask1 = (qqqE > 0) & ~np.isnan(thetab)
plt.figure(figsize=(7,6))
plt.hist2d(thetab[mask1], qqqE[mask1], bins=200)
plt.xlabel("thetab")
plt.ylabel("E-QQQ")
plt.title(f"{particle} ({treename}) QQQ Energy vs Theta")
plt.colorbar(label="Counts")
#plt.xlim(0,30)
#plt.ylim(0,0.45)
plt.tight_layout()
plt.savefig(f"{base}/QQQE_vs_theta.png", dpi=300)
plt.show()
plt.figure(figsize=(7,5)) plt.figure(figsize=(7,5))
plt.hist(Elost, bins=200) plt.hist(Elost, bins=200)
@ -1014,18 +1139,6 @@ class MyInteractiveApp(cmd.Cmd):
plt.legend(loc='upper right') plt.legend(loc='upper right')
plt.show() 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")
plt.figure(figsize=(7,5)) plt.figure(figsize=(7,5))
plt.hist(dA, bins=100, label='dA', histtype='step') plt.hist(dA, bins=100, label='dA', histtype='step')
@ -1053,8 +1166,9 @@ class MyInteractiveApp(cmd.Cmd):
print("Value error") print("Value error")
try: try:
mask = Esx3 > 0
plt.figure(figsize=(7,6)) plt.figure(figsize=(7,6))
plt.hist2d(EA, Esx3, bins=200) plt.hist2d(EA[mask], Esx3[mask], bins=200)
plt.xlabel("EA (MeV)") plt.xlabel("EA (MeV)")
plt.ylabel("Esx3 (MeV)") plt.ylabel("Esx3 (MeV)")
plt.title(f"{particle} ({treename}) Anode vs SX3 Energy") plt.title(f"{particle} ({treename}) Anode vs SX3 Energy")
@ -1081,7 +1195,7 @@ class MyInteractiveApp(cmd.Cmd):
try: try:
mask1 = ~np.isnan(Esx3) & (Esx3 > 0) mask1 = ~np.isnan(Esx3) & (Esx3 > 0)
plt.figure(figsize=(7,6)) plt.figure(figsize=(7,6))
plt.hist2d(Esx3[mask1], Eprop[mask1], bins=200) plt.hist2d(Esx3[mask1], Eprop[mask1] * data["thetab"][mask1], bins=200)
plt.ylabel("PCEnergy") plt.ylabel("PCEnergy")
plt.xlabel("SX3 Energy (MeV)") plt.xlabel("SX3 Energy (MeV)")
plt.title(f"{particle} ({treename}) Energy Propagation Difference vs SX3 Energy") plt.title(f"{particle} ({treename}) Energy Propagation Difference vs SX3 Energy")
@ -1094,21 +1208,10 @@ class MyInteractiveApp(cmd.Cmd):
except ValueError: except ValueError:
print("Value error") print("Value error")
try: #try:
mask1 = ~np.isnan(Ei) & (Esx3 > 0)
plt.figure(figsize=(7,6)) #except ValueError:
plt.hist2d(data["thetab"], Ei, bins=200) # print("Value error")
plt.xlabel("thetab")
plt.ylabel("Event Energy (MeV)")
plt.title(f"{particle} ({treename}) Energy vs Theta")
plt.colorbar(label="Counts")
#plt.xlim(0,30)
#plt.ylim(0,0.45)
plt.tight_layout()
plt.savefig(f"{base}/E_vs_theta.png", dpi=300)
plt.show()
except ValueError:
print("Value error")
branch_names = [] branch_names = []
@ -1126,9 +1229,8 @@ class MyInteractiveApp(cmd.Cmd):
all_branches = self.tree.arrays(branch_names, library="np", entry_stop=max_events) all_branches = self.tree.arrays(branch_names, library="np", entry_stop=max_events)
# Keep only events with thetab >= 45 # Keep only events with thetab >= 45
#thetab_mask = all_branches["thetab"] >= 0 #mask = (all_branches["Tb"] > 0) & (all_branches["qqqTb"] > 0)
for branch in branch_names: for branch in branch_names:
values = all_branches[branch] values = all_branches[branch]
@ -1143,7 +1245,7 @@ class MyInteractiveApp(cmd.Cmd):
continue continue
# Apply the thetab cut # Apply the thetab cut
values = values #values = values[mask]
values = values[~np.isnan(values)] values = values[~np.isnan(values)]
@ -1164,6 +1266,8 @@ class MyInteractiveApp(cmd.Cmd):
plt.close() plt.close()
else: else:
print("No branches found to histogram.") print("No branches found to histogram.")
print(f"Plotting complete ({treename}).") print(f"Plotting complete ({treename}).")
@ -1408,8 +1512,8 @@ class MyInteractiveApp(cmd.Cmd):
plt.margins(0) plt.margins(0)
plt.xlabel("SX3 Energy (MeV)") plt.xlabel("SX3 Energy (MeV)")
plt.ylabel("PCEnergy x Sin(theta)") plt.ylabel("PCEnergy x Sin(theta)")
plt.xlim(0,35) #plt.xlim(0,35)
plt.ylim(0,0.6) #plt.ylim(0,0.6)
cbar = plt.colorbar() cbar = plt.colorbar()
cbar.set_label("Counts") cbar.set_label("Counts")
plt.minorticks_on() plt.minorticks_on()