sausage eliminated
This commit is contained in:
James Szalkie
parent
d94795ae33
commit
258306d442
|
|
@ -57,8 +57,8 @@ private:
|
|||
const int numDet = 12;
|
||||
const float radius = 88;
|
||||
const float width = 40;
|
||||
const float length = 75;
|
||||
const float gap = 46;
|
||||
const float length = 174.3; // 75
|
||||
const float gap = 0; // 46
|
||||
|
||||
short id; // -1 when no hit
|
||||
short chUp;
|
||||
|
|
|
|||
|
|
@ -96,13 +96,15 @@ int main(int argc, char **argv){
|
|||
if( argc >= 2 ) numEvent = atoi(argv[1]);
|
||||
TransferReaction transfer;
|
||||
|
||||
transfer.SetA(14, 7, 0); // 18Ne projectile
|
||||
transfer.SetIncidentEnergyAngle((42.82/14), 0, 0); // KEA in MeV/u, theta and phi in rad
|
||||
transfer.SetA(18, 10, 0); // 18Ne projectile
|
||||
//transfer.SetIncidentEnergyAngle(0, 0, 0); // KEA in MeV/u, theta and phi in rad
|
||||
TGraph* elossBeam = LoadELoss("../ELoss/E_vs_x_Ne-18.dat");
|
||||
transfer.Seta(4, 2); // 4He target
|
||||
transfer.Setb(1, 1); // outgoing proton from the primary transfer
|
||||
transfer.SetB(17, 8); // 21Na* heavy product
|
||||
transfer.SetB(21, 11); // 21Na* heavy product
|
||||
const double beamA = 18.0; // mass number of 18Ne beam
|
||||
|
||||
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;
|
||||
|
|
@ -110,12 +112,14 @@ int main(int argc, char **argv){
|
|||
|
||||
// Excited state lists (projectile and heavy-product excitation states)
|
||||
std::vector<float> ExAList = {0}; // 18Ne projectile excitations in MeV
|
||||
std::vector<float> ExList = {0}; // 21Na* excitation in MeV
|
||||
std::vector<float> ExList = {2.797, 2.829, 3.544, 3.862, 4.294}; // 21Na* excitation 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};
|
||||
double vertexXRange[2] = { -5, 5}; // mm - 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)
|
||||
const double beamEntranceZ = -174.3; //vertexZRange[0]; // mm, assumed beam entrance into the gas
|
||||
|
||||
|
||||
// detector resolution / uncertainty parameters
|
||||
double sigmaSX3_W = 0; // mm, if < 0 use mid-point (no spread in SX3 horizontal dimension)
|
||||
|
|
@ -168,8 +172,17 @@ int main(int argc, char **argv){
|
|||
// beam and CM variables saved in tree
|
||||
double KEA;
|
||||
double KEA2;
|
||||
double beamPath_cm;
|
||||
double beamEnergy;
|
||||
double beamEnergyLoss;
|
||||
tree1->Branch("beamKEA", &KEA, "beamKEA/D");
|
||||
tree2->Branch("beamKEA", &KEA2, "beamKEA/D");
|
||||
tree1->Branch("beamPath_cm", &beamPath_cm, "beamPath_cm/D");
|
||||
tree2->Branch("beamPath_cm", &beamPath_cm, "beamPath_cm/D");
|
||||
tree1->Branch("beamEnergy", &beamEnergy, "beamEnergy/D");
|
||||
tree2->Branch("beamEnergy", &beamEnergy, "beamEnergy/D");
|
||||
tree1->Branch("beamEnergyLoss", &beamEnergyLoss, "beamEnergyLoss/D");
|
||||
tree2->Branch("beamEnergyLoss", &beamEnergyLoss, "beamEnergyLoss/D");
|
||||
|
||||
double thetaCM, phiCM;
|
||||
double thetaCM2, phiCM2;
|
||||
|
|
@ -319,6 +332,22 @@ int main(int argc, char **argv){
|
|||
// recalc kinematic constants for chosen states
|
||||
transfer.CalReactionConstant();
|
||||
|
||||
// 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);
|
||||
|
||||
// compute beam energy at the event vertex from the gas path length
|
||||
beamPath_cm = TVector3(vertexZ - beamEntranceZ, vertexX, vertexY).Mag() * 0.1;
|
||||
if( beamPath_cm < 0 ) beamPath_cm = 0;
|
||||
beamEnergy = elossBeam->Eval(beamPath_cm); // MeV
|
||||
beamEnergyLoss = elossBeam->Eval(0.0) - beamEnergy;
|
||||
KEA = beamEnergy / beamA;
|
||||
transfer.SetIncidentEnergyAngle(KEA, 0, 0);
|
||||
transfer.CalReactionConstant();
|
||||
|
||||
// isotropic CM direction
|
||||
thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ; // polar angle in CM frame
|
||||
phiCM = (gRandom->Rndm() - 0.5) * TMath::TwoPi();
|
||||
|
|
@ -364,13 +393,6 @@ int main(int argc, char **argv){
|
|||
|
||||
delete [] output;
|
||||
|
||||
// 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());
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
501
ELoss/E_vs_x_Ne-18.dat
Normal file
501
ELoss/E_vs_x_Ne-18.dat
Normal file
|
|
@ -0,0 +1,501 @@
|
|||
Distance_cm Energy_MeV
|
||||
-0.0 80.0
|
||||
0.16713525463015264 79.83969939879759
|
||||
0.33406473805004305 79.6793987975952
|
||||
0.5007884444544184 79.51909819639279
|
||||
0.6673063680316278 79.35879759519038
|
||||
0.8336185029635974 79.19849699398797
|
||||
0.9997248434257149 79.03819639278558
|
||||
1.1656253835868005 78.87789579158317
|
||||
1.3313201176089287 78.71759519038076
|
||||
1.4968090396473945 78.55729458917835
|
||||
1.662092143850605 78.39699398797595
|
||||
1.8271694243599683 78.23669338677355
|
||||
1.9920408753098522 78.07639278557114
|
||||
2.1567064908273967 77.91609218436874
|
||||
2.3211662650324683 77.75579158316633
|
||||
2.485420192037525 77.59549098196393
|
||||
2.6494682659475672 77.43519038076153
|
||||
2.813310480859938 77.27488977955912
|
||||
2.97694683086427 77.11458917835671
|
||||
3.1403773100423384 76.95428857715432
|
||||
3.303601912468001 76.79398797595191
|
||||
3.466620632206991 76.6336873747495
|
||||
3.6294334633168486 76.47338677354709
|
||||
3.7920403998467678 76.3130861723447
|
||||
3.9544414358375213 76.15278557114229
|
||||
4.116636565321244 75.99248496993988
|
||||
4.278625782321351 75.83218436873747
|
||||
4.440409080852373 75.67188376753508
|
||||
4.601986454919871 75.51158316633267
|
||||
4.763357898520203 75.35128256513026
|
||||
4.924523405640434 75.19098196392785
|
||||
5.085482970258169 75.03068136272545
|
||||
5.246236586341382 74.87038076152305
|
||||
5.406784247848318 74.71008016032064
|
||||
5.567125948727237 74.54977955911824
|
||||
5.727261682916312 74.38947895791583
|
||||
5.8871914443434274 74.22917835671343
|
||||
6.0469152269260595 74.06887775551102
|
||||
6.206433024571017 73.90857715430862
|
||||
6.365744831174311 73.74827655310621
|
||||
6.524850640620946 73.58797595190381
|
||||
6.683750446784783 73.4276753507014
|
||||
6.842444243528258 73.267374749499
|
||||
7.0009320247022435 73.10707414829659
|
||||
7.1592137841458126 72.9467735470942
|
||||
7.31728951568609 72.78647294589179
|
||||
7.475159213137954 72.62617234468938
|
||||
7.632822870303875 72.46587174348697
|
||||
7.790280480973665 72.30557114228458
|
||||
7.947532038924305 72.14527054108217
|
||||
8.10457753791963 71.98496993987976
|
||||
8.261416971710144 71.82466933867735
|
||||
8.418050334032769 71.66436873747494
|
||||
8.574477618610576 71.50406813627255
|
||||
8.73069881915259 71.34376753507014
|
||||
8.886713929353444 71.18346693386773
|
||||
9.042522942893175 71.02316633266533
|
||||
9.198125853436917 70.86286573146293
|
||||
9.353522654634682 70.70256513026052
|
||||
9.508713340120996 70.54226452905812
|
||||
9.66369790351466 70.38196392785571
|
||||
9.818476338418428 70.22166332665331
|
||||
9.973048638418758 70.0613627254509
|
||||
10.127414797085422 69.9010621242485
|
||||
10.28157480797124 69.74076152304609
|
||||
10.43552866461173 69.5804609218437
|
||||
10.589276360524838 69.42016032064129
|
||||
10.742817889210507 69.25985971943888
|
||||
10.896153244150392 69.09955911823647
|
||||
11.049282418807483 68.93925851703408
|
||||
11.202205406625788 68.77895791583167
|
||||
11.354922201029893 68.61865731462926
|
||||
11.50743279542463 68.45835671342685
|
||||
11.659737183194682 68.29805611222444
|
||||
11.811835357704176 68.13775551102205
|
||||
11.96372731229633 67.97745490981964
|
||||
12.115413040292964 67.81715430861723
|
||||
12.266892534994133 67.65685370741483
|
||||
12.418165789677664 67.49655310621243
|
||||
12.569232797598776 67.33625250501002
|
||||
12.72009355198955 67.17595190380761
|
||||
12.870748046058516 67.0156513026052
|
||||
13.021196272990165 66.85535070140281
|
||||
13.17143822594452 66.6950501002004
|
||||
13.321473898056563 66.534749498998
|
||||
13.471303282435787 66.37444889779559
|
||||
13.62092637216566 66.2141482965932
|
||||
13.770343160303154 66.05384769539079
|
||||
13.919553639878124 65.89354709418838
|
||||
14.068557803892825 65.73324649298597
|
||||
14.217355645321316 65.57294589178358
|
||||
14.365947157108948 65.41264529058117
|
||||
14.514332332171698 65.25234468937876
|
||||
14.662511163395624 65.09204408817635
|
||||
14.810483643636248 64.93174348697394
|
||||
14.958249765717909 64.77144288577155
|
||||
15.105809522433189 64.61114228456914
|
||||
15.253162906542174 64.45084168336673
|
||||
15.400309910771849 64.29054108216432
|
||||
15.547250527815386 64.13024048096193
|
||||
15.69398475033151 63.969939879759515
|
||||
15.840512570943682 63.80963927855711
|
||||
15.986833982239466 63.649338677354706
|
||||
16.132948976769757 63.4890380761523
|
||||
16.278857547048005 63.3287374749499
|
||||
16.42455968554947 63.16843687374749
|
||||
16.570055384710404 63.00813627254509
|
||||
16.71534463692726 62.84783567134268
|
||||
16.86042743455584 62.68753507014028
|
||||
17.005303769910494 62.52723446893787
|
||||
17.149973635263198 62.36693386773547
|
||||
17.294437022842725 62.20663326653306
|
||||
17.438693924833697 62.04633266533066
|
||||
17.5827443333757 61.88603206412825
|
||||
17.726588240562307 61.72573146292585
|
||||
17.870225638440154 61.56543086172344
|
||||
18.013656519007906 61.405130260521034
|
||||
18.156880874215283 61.24482965931863
|
||||
18.299898695962035 61.084529058116225
|
||||
18.442709976096864 60.92422845691382
|
||||
18.58531470641638 60.763927855711415
|
||||
18.72771287866396 60.603627254509014
|
||||
18.869904484528682 60.443326653306606
|
||||
19.01188951564411 60.283026052104205
|
||||
19.15366796358719 60.1227254509018
|
||||
19.295239819876986 59.962424849699396
|
||||
19.436605075973514 59.80212424849699
|
||||
19.57776372327643 59.641823647294586
|
||||
19.718715753123806 59.48152304609218
|
||||
19.859461156790765 59.32122244488978
|
||||
19.999999925488197 59.16092184368737
|
||||
20.140332050361344 59.00062124248497
|
||||
20.28045752248844 58.84032064128256
|
||||
20.420376332879247 58.68002004008016
|
||||
20.560088472473627 58.51971943887775
|
||||
20.69959393214001 58.35941883767535
|
||||
20.83889270267392 58.19911823647294
|
||||
20.977984774796354 58.03881763527053
|
||||
21.116870139152223 57.87851703406813
|
||||
21.25554878630873 57.718216432865724
|
||||
21.394020706753658 57.55791583166332
|
||||
21.53228589089372 57.397615230460914
|
||||
21.670344329052767 57.23731462925851
|
||||
21.808196011470045 57.077014028056105
|
||||
21.945840928298335 56.916713426853704
|
||||
22.083279069602135 56.756412825651296
|
||||
22.220510425355712 56.596112224448895
|
||||
22.357534985441198 56.43581162324649
|
||||
22.49435273964655 56.275511022044086
|
||||
22.630963677663576 56.11521042084168
|
||||
22.767367789085785 55.954909819639276
|
||||
22.90356506340633 55.79460921843687
|
||||
23.03955549001576 55.63430861723447
|
||||
23.175339058199874 55.47400801603206
|
||||
23.310915757137373 55.31370741482966
|
||||
23.446285575897594 55.15340681362725
|
||||
23.58144850343808 54.99310621242485
|
||||
23.716404528602197 54.83280561122244
|
||||
23.851153640116614 54.67250501002003
|
||||
23.985695826588763 54.51220440881763
|
||||
24.12003107650427 54.35190380761522
|
||||
24.254159378224237 54.19160320641282
|
||||
24.388080719982597 54.03130260521041
|
||||
24.52179508988326 53.87100200400801
|
||||
24.655302475897322 53.710701402805604
|
||||
24.788602865860124 53.5504008016032
|
||||
24.921696247468297 53.390100200400795
|
||||
25.0545826082767 53.229799599198394
|
||||
25.187261935695332 53.069498997995986
|
||||
25.319734216986117 52.909198396793585
|
||||
25.45199943925968 52.748897795591176
|
||||
25.58405758947199 52.588597194388775
|
||||
25.715908654420975 52.42829659318637
|
||||
25.847552620743016 52.267995991983966
|
||||
25.978989474909405 52.10769539078156
|
||||
26.11021920322267 51.94739478957916
|
||||
26.24124179181288 51.78709418837675
|
||||
26.372057226633782 51.62679358717435
|
||||
26.502665493458952 51.46649298597194
|
||||
26.633066577877752 51.30619238476954
|
||||
26.763260465291278 51.14589178356713
|
||||
26.893247140908148 50.98559118236472
|
||||
27.023026589740244 50.82529058116232
|
||||
27.152598796598348 50.66498997995991
|
||||
27.281963746087623 50.50468937875751
|
||||
27.411121422603074 50.3443887775551
|
||||
27.540071810324818 50.1840881763527
|
||||
27.668814893213327 50.023787575150294
|
||||
27.797350655004475 49.86348697394789
|
||||
27.92567907920456 49.703186372745485
|
||||
28.05380014908511 49.542885771543084
|
||||
28.18171384767767 49.382585170340676
|
||||
28.309420157768376 49.222284569138274
|
||||
28.43691906189247 49.061983967935866
|
||||
28.564210542328635 48.901683366733465
|
||||
28.69129458109325 48.74138276553106
|
||||
28.818171159934437 48.581082164328656
|
||||
28.944840260326067 48.42078156312625
|
||||
29.071301863461503 48.26048096192385
|
||||
29.197555950247317 48.10018036072144
|
||||
29.323602501296758 47.93987975951904
|
||||
29.449441496923143 47.77957915831663
|
||||
29.57507291713302 47.61927855711422
|
||||
29.700496741619226 47.45897795591182
|
||||
29.82571294975377 47.29867735470941
|
||||
29.95072152058049 47.13837675350701
|
||||
30.075522432807634 46.9780761523046
|
||||
30.20011566480015 46.8177755511022
|
||||
30.32450119457191 46.65747494989979
|
||||
30.448678999777613 46.49717434869739
|
||||
30.57264905770465 46.336873747494984
|
||||
30.696411345264607 46.17657314629258
|
||||
30.819965838984718 46.016272545090175
|
||||
30.94331251499898 45.85597194388777
|
||||
31.066451349039156 45.695671342685365
|
||||
31.189382316425476 45.535370741482964
|
||||
31.3121053920572 45.375070140280556
|
||||
31.43462055040285 45.214769539078155
|
||||
31.556927765490304 45.05446893787575
|
||||
31.679027010896565 44.894168336673346
|
||||
31.800918259737347 44.73386773547094
|
||||
31.92260148465635 44.57356713426854
|
||||
32.04407665781432 44.41326653306613
|
||||
32.16534375087781 44.25296593186372
|
||||
32.28640273500765 44.09266533066132
|
||||
32.4072535808472 43.93236472945891
|
||||
32.52789625851022 43.77206412825651
|
||||
32.64833073756853 43.6117635270541
|
||||
32.76855698703928 43.4514629258517
|
||||
32.888574975372 43.29116232464929
|
||||
33.00838467043521 43.13086172344689
|
||||
33.12798603950282 42.97056112224448
|
||||
33.247379049240074 42.81026052104208
|
||||
33.36656366568924 42.649959919839674
|
||||
33.48553985425489 42.48965931863727
|
||||
33.60430757968881 42.329358717434864
|
||||
33.72286680607454 42.16905811623246
|
||||
33.841217496811566 42.008757515030055
|
||||
33.959359614598995 41.848456913827654
|
||||
34.07729312141897 41.688156312625246
|
||||
34.19501797851952 41.527855711422845
|
||||
34.31253414639713 41.36755511022044
|
||||
34.42984158477867 41.207254509018036
|
||||
34.54694025260308 41.04695390781563
|
||||
34.66383010800239 40.88665330661322
|
||||
34.78051110828239 40.72635270541082
|
||||
34.896983209902785 40.56605210420841
|
||||
35.013246368456734 40.40575150300601
|
||||
35.12930053865002 40.2454509018036
|
||||
35.245145674279534 40.0851503006012
|
||||
35.360781728211336 39.92484969939879
|
||||
35.47620865235801 39.76454909819639
|
||||
35.59142639765558 39.60424849699398
|
||||
35.706434914039676 39.44394789579158
|
||||
35.82123415042119 39.28364729458917
|
||||
35.93582405466124 39.12334669338677
|
||||
36.050204573545514 38.96304609218436
|
||||
36.16437565275787 38.80274549098196
|
||||
36.27833723685334 38.642444889779554
|
||||
36.39208926923033 38.48214428857715
|
||||
36.50563169210217 38.321843687374745
|
||||
36.61896444646781 38.161543086172344
|
||||
36.73208747208184 38.001242484969936
|
||||
36.84500070742364 37.840941883767535
|
||||
36.95770408966574 37.68064128256513
|
||||
37.070197554641325 37.52034068136272
|
||||
37.18248103681086 37.36004008016032
|
||||
37.294554469227855 37.19973947895791
|
||||
37.40641778350363 37.03943887775551
|
||||
37.51807090977123 36.8791382765531
|
||||
37.62951377664829 36.7188376753507
|
||||
37.740746311198905 36.55853707414829
|
||||
37.851768438894496 36.39823647294589
|
||||
37.96258008357358 36.23793587174348
|
||||
38.07318116740042 36.07763527054108
|
||||
38.18357161082264 35.91733466933867
|
||||
38.29375133252751 35.75703406813627
|
||||
38.40372024939723 35.59673346693386
|
||||
38.51347827646278 35.43643286573146
|
||||
38.623025326856656 35.27613226452905
|
||||
38.732361311764194 35.11583166332665
|
||||
38.841486140373604 34.955531062124244
|
||||
38.950399719824546 34.79523046092184
|
||||
39.05910195515535 34.634929859719435
|
||||
39.16759274924867 34.474629258517034
|
||||
39.275872002775714 34.314328657314626
|
||||
39.3839396141388 34.15402805611222
|
||||
39.491795479412346 33.993727454909816
|
||||
39.59943949228223 33.83342685370741
|
||||
39.70687154398332 33.67312625250501
|
||||
39.81409152323534 33.5128256513026
|
||||
39.92109931617683 33.3525250501002
|
||||
40.02789480629727 33.19222444889779
|
||||
40.13447787436722 33.03192384769539
|
||||
40.24084839836649 32.87162324649298
|
||||
40.34700625341019 32.71132264529058
|
||||
40.45295131167273 32.55102204408817
|
||||
40.55868344230947 32.39072144288577
|
||||
40.6642025113763 32.23042084168336
|
||||
40.76950838174669 32.07012024048096
|
||||
40.87460091302642 31.909819639278556
|
||||
40.97947996146579 31.74951903807615
|
||||
41.08414537986927 31.589218436873747
|
||||
41.18859701750243 31.428917835671342
|
||||
41.292834719996215 31.268617234468937
|
||||
41.396858329248325 31.108316633266533
|
||||
41.50066768332168 30.948016032064128
|
||||
41.604262616339874 30.787715430861724
|
||||
41.707642958379516 30.62741482965932
|
||||
41.81080853535932 30.467114228456914
|
||||
41.91375916892593 30.30681362725451
|
||||
42.01649467633623 30.146513026052105
|
||||
42.11901487033617 29.9862124248497
|
||||
42.22131955903592 29.825911823647296
|
||||
42.32340854578118 29.66561122244489
|
||||
42.42528162902067 29.505310621242486
|
||||
42.52693860216951 29.345010020040082
|
||||
42.628379253468445 29.184709418837677
|
||||
42.72960336583878 29.02440881763527
|
||||
42.830610716732785 28.864108216432864
|
||||
42.9314010779796 28.70380761523046
|
||||
43.03197421562627 28.543507014028055
|
||||
43.13232988977391 28.38320641282565
|
||||
43.23246785440877 28.222905811623246
|
||||
43.33238785722805 28.06260521042084
|
||||
43.43208963946024 27.902304609218437
|
||||
43.53157293567985 27.742004008016032
|
||||
43.63083747361634 27.581703406813627
|
||||
43.72988297395699 27.421402805611223
|
||||
43.828709150143595 27.261102204408818
|
||||
43.927315708162666 27.100801603206413
|
||||
44.02570234632903 26.94050100200401
|
||||
44.123868755062475 26.780200400801604
|
||||
44.221814616657355 26.6198997995992
|
||||
44.31953960504472 26.459599198396795
|
||||
44.4170433855469 26.29929859719439
|
||||
44.514325614624134 26.138997995991986
|
||||
44.611385939613086 25.97869739478958
|
||||
44.70822399845685 25.818396793587176
|
||||
44.804839419426244 25.65809619238477
|
||||
44.901231820832024 25.497795591182363
|
||||
44.99740081072773 25.33749498997996
|
||||
45.09334598660287 25.177194388777554
|
||||
45.189066935065966 25.01689378757515
|
||||
45.28456323151731 24.856593186372745
|
||||
45.37983443981088 24.69629258517034
|
||||
45.47488011190513 24.535991983967936
|
||||
45.56969978750229 24.37569138276553
|
||||
45.66429299367563 24.215390781563126
|
||||
45.75865924448445 24.05509018036072
|
||||
45.85279804057618 23.894789579158317
|
||||
45.94670886877532 23.734488977955913
|
||||
46.04039120165851 23.574188376753508
|
||||
46.13384449711553 23.413887775551103
|
||||
46.227068197895484 23.2535871743487
|
||||
46.320061731137756 23.093286573146294
|
||||
46.41282450788721 22.93298597194389
|
||||
46.50535592259301 22.772685370741485
|
||||
46.59765535259054 22.61238476953908
|
||||
46.68972828343022 22.452084168336675
|
||||
46.781591894934756 22.29178356713427
|
||||
46.873258589408394 22.131482965931863
|
||||
46.964729457462816 21.971182364729458
|
||||
47.05600463342356 21.810881763527053
|
||||
47.14708425265793 21.65058116232465
|
||||
47.2379684515823 21.490280561122244
|
||||
47.32865736767084 21.32997995991984
|
||||
47.41915113946602 21.169679358717435
|
||||
47.509449906590845 21.00937875751503
|
||||
47.59955380976344 20.849078156312626
|
||||
47.689462990813915 20.68877755511022
|
||||
47.779177592704045 20.528476953907816
|
||||
47.86869775955004 20.36817635270541
|
||||
47.95802363664879 20.207875751503007
|
||||
48.04715537050797 20.047575150300602
|
||||
48.13609310888054 19.887274549098198
|
||||
48.22483700080402 19.726973947895793
|
||||
48.31338719664531 19.56667334669339
|
||||
48.4017438481514 19.406372745490984
|
||||
48.48990710850689 19.24607214428858
|
||||
48.577877132398925 19.085771543086175
|
||||
48.665654076090405 18.92547094188377
|
||||
48.7532380975024 18.76517034068136
|
||||
48.84062935630667 18.604869739478957
|
||||
48.927828014029565 18.444569138276552
|
||||
49.01483423416827 18.284268537074148
|
||||
49.10164818232089 18.123967935871743
|
||||
49.18827002633179 17.96366733466934
|
||||
49.27469993645373 17.803366733466934
|
||||
49.36093808552857 17.64306613226453
|
||||
49.446984649188565 17.482765531062125
|
||||
49.53283980608019 17.32246492985972
|
||||
49.618503738113006 17.162164328657315
|
||||
49.70397663073599 17.00186372745491
|
||||
49.78925867324422 16.841563126252506
|
||||
49.87435005911901 16.6812625250501
|
||||
49.95925098640481 16.520961923847697
|
||||
50.043961658126726 16.360661322645292
|
||||
50.12848228275268 16.200360721442888
|
||||
50.212813074704805 16.040060120240483
|
||||
50.29695425492496 15.879759519038075
|
||||
50.380906051499935 15.71945891783567
|
||||
50.46466870035226 15.559158316633265
|
||||
50.54824244600334 15.39885771543086
|
||||
50.63162754241612 15.238557114228456
|
||||
50.714824253925336 15.078256513026052
|
||||
50.797832856264215 14.917955911823647
|
||||
50.88065363769728 14.757655310621242
|
||||
50.96328690027001 14.597354709418838
|
||||
51.04573296118719 14.437054108216431
|
||||
51.127992154332844 14.276753507014027
|
||||
51.21006483194625 14.116452905811622
|
||||
51.29195136646969 13.956152304609217
|
||||
51.373652152585464 13.795851703406813
|
||||
51.45516760946135 13.635551102204408
|
||||
51.53649818322577 13.475250501002003
|
||||
51.617644349696086 13.314949899799599
|
||||
51.69860661738592 13.154649298597194
|
||||
51.77938553082005 12.99434869739479
|
||||
51.859997377378576 12.834048096192385
|
||||
51.940470139903056 12.673747494989978
|
||||
52.020817710343564 12.513446893787574
|
||||
52.101043131727295 12.35314629258517
|
||||
52.18114871928441 12.192845691382765
|
||||
52.26113691037399 12.03254509018036
|
||||
52.341010274159125 11.872244488977955
|
||||
52.42077152218919 11.71194388777555
|
||||
52.50042351998291 11.551643286573146
|
||||
52.57996929971611 11.391342685370741
|
||||
52.659412074129975 11.231042084168337
|
||||
52.73875525178858 11.07074148296593
|
||||
52.8180024538295 10.910440881763526
|
||||
52.89715753236777 10.750140280561121
|
||||
52.9762245907325 10.589839679358716
|
||||
53.05520800573647 10.429539078156312
|
||||
53.1341124522032 10.269238476953907
|
||||
53.21294293000315 10.108937875751502
|
||||
53.291704793881614 9.948637274549098
|
||||
53.37040378639588 9.788336673346693
|
||||
53.44904607431941 9.628036072144289
|
||||
53.527638288916116 9.467735470941884
|
||||
53.60618757054037 9.307434869739478
|
||||
53.68470161807782 9.147134268537073
|
||||
53.76318874381103 8.986833667334668
|
||||
53.841657934373046 8.826533066132264
|
||||
53.920118918543174 8.666232464929859
|
||||
53.99858224274524 8.505931863727454
|
||||
54.07705935523119 8.34563126252505
|
||||
54.1555627000758 8.185330661322645
|
||||
54.23410582227539 8.02503006012024
|
||||
54.31270348543876 7.864729458917835
|
||||
54.39137180378868 7.70442885771543
|
||||
54.47012839046343 7.544128256513026
|
||||
54.54899252442924 7.383827655310621
|
||||
54.62798533869604 7.2235270541082155
|
||||
54.70713003298487 7.063226452905811
|
||||
54.78639077236708 6.902925851703406
|
||||
54.8656671752427 6.742625250501002
|
||||
54.94490987977844 6.582324649298597
|
||||
55.02413071285859 6.422024048096192
|
||||
55.10334829033552 6.261723446893787
|
||||
55.18258344018353 6.101422845691382
|
||||
55.26185947485493 5.9411222444889775
|
||||
55.34120250050544 5.780821643286573
|
||||
55.420641768914145 5.620521042084168
|
||||
55.50021007903053 5.460220440881763
|
||||
55.579944236439935 5.299919839679358
|
||||
55.65988558071507 5.1396192384769535
|
||||
55.74008059270073 4.979318637274549
|
||||
55.82058159637453 4.819018036072144
|
||||
55.90144757318593 4.658717434869739
|
||||
55.98274511089628 4.498416833667334
|
||||
56.06454951418689 4.338116232464929
|
||||
56.14694611102537 4.177815631262525
|
||||
56.23003179746609 4.01751503006012
|
||||
56.313916874873435 3.857214428857715
|
||||
56.39872724841941 3.6969138276553104
|
||||
56.48460707541674 3.5366132264529053
|
||||
56.57172197844761 3.3763126252505007
|
||||
56.66026297398809 3.216012024048096
|
||||
56.75045131617868 3.055711422845691
|
||||
56.842544523282776 2.8954108216432863
|
||||
56.936843949828585 2.7351102204408813
|
||||
57.03370440364967 2.5748096192384766
|
||||
57.133546504624924 2.414509018036072
|
||||
57.23687277369295 2.254208416833667
|
||||
57.344288880173586 2.0939078156312623
|
||||
57.4565321519782 1.9336072144288576
|
||||
57.57451052057477 1.7733066132264528
|
||||
57.69935680289275 1.6130060120240481
|
||||
57.832506115155404 1.4527054108216433
|
||||
57.975809220578434 1.2924048096192384
|
||||
58.13170362498847 1.1321042084168336
|
||||
58.30348118145025 0.9718036072144288
|
||||
58.49572437647055 0.8115030060120241
|
||||
58.715052349153865 0.6512024048096192
|
||||
58.97145636708799 0.4909018036072144
|
||||
59.276292273963826 0.3306012024048096
|
||||
59.64908939725906 0.1703006012024048
|
||||
60.08468062361392 0.01
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
|
|
@ -20,7 +20,6 @@ import os
|
|||
import periodictable as pt
|
||||
import re
|
||||
from matplotlib.colors import LinearSegmentedColormap
|
||||
import matplotlib as mpl
|
||||
|
||||
# ROOT-like styling
|
||||
plt.rcParams.update({
|
||||
|
|
@ -79,9 +78,8 @@ interp_cache = {}
|
|||
def make_E_vs_x(z, mass_u, emax_mev, label, npoints, P_TORR, TEMP_K):
|
||||
# GAS SETUP
|
||||
R = 8.3144
|
||||
P_TORR = 400
|
||||
TEMP_K = 293.15
|
||||
R = 8.3144
|
||||
P_TORR = P_TORR
|
||||
TEMP_K = TEMP_K
|
||||
|
||||
# Gas density
|
||||
p_pa = P_TORR * 133.322
|
||||
|
|
@ -222,8 +220,6 @@ def energy_loss(particle, Ei, dl):
|
|||
|
||||
return np.maximum(Ef, 0.0)
|
||||
|
||||
#def energy_loss_angular(particle, Ei, theta):
|
||||
|
||||
def energy_reconstruction(particle, Ef, dl):
|
||||
|
||||
E_of_x, x_of_E = get_interpolators(particle)
|
||||
|
|
@ -242,11 +238,10 @@ def energy_distance(particle, Ei, Ef):
|
|||
|
||||
xi = x_of_E(Ei)
|
||||
|
||||
xf = x_of_E(Ef)
|
||||
xf = x_of_E(Ef) * .686
|
||||
|
||||
return np.abs(xf - xi)
|
||||
|
||||
|
||||
def resolve_particle(name):
|
||||
name = name.lower().strip().rstrip("s")
|
||||
if name in particles:
|
||||
|
|
@ -274,7 +269,6 @@ def resolve_particle(name):
|
|||
except Exception:
|
||||
raise ValueError(f"Unknown particle/isotope: {name}")
|
||||
|
||||
|
||||
def update_plot_data(name, values):
|
||||
for i, (existing_name, _) in enumerate(plot_data):
|
||||
if existing_name == name:
|
||||
|
|
@ -296,19 +290,17 @@ def calculate_distance_tree2(vz, theta, z_max=34.86):
|
|||
cos_theta = np.where(np.abs(cos_theta) < 1e-10, 1e-10, cos_theta)
|
||||
return np.abs((z_max - vz) / cos_theta)
|
||||
|
||||
|
||||
def load_tree_arrays(tree, treename, max_events=None):
|
||||
branches = ["Tb", "thetab", "sx3Z", "vX", "vY", "vZ", "sx3ID", "aID", "cID", "aDist", "cDist"]
|
||||
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)
|
||||
|
||||
mask = data["thetab"] >= 0
|
||||
data = {key: value[mask] for key, value in data.items()}
|
||||
#mask = data["thetab"] >= 0
|
||||
#data = {key: value[mask] for key, value in data.items()}
|
||||
Ei = data["Tb"]
|
||||
theta = np.radians(data["thetab"])
|
||||
vX = data["vX"]
|
||||
|
|
@ -340,7 +332,7 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
|
|||
|
||||
sin_theta = np.sin(theta)
|
||||
sin_theta = np.where(sin_theta != 0, sin_theta, 1e-10)
|
||||
radii = np.array([3.7, 4.2])
|
||||
radii = np.array([3.7, 4.3])
|
||||
dA = (radii[0] - np.sqrt((vX/10)**2 + (vY/10)**2))/ sin_theta
|
||||
dC = (radii[1] - np.sqrt((vX/10)**2 + (vY/10)**2))/ sin_theta
|
||||
|
||||
|
|
@ -392,10 +384,9 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
|
|||
"Eprop": Eprop,
|
||||
"dA": dA,
|
||||
"dC": dC,
|
||||
"thetab": data["thetab"]
|
||||
"thetab": np.degrees(theta)
|
||||
}
|
||||
|
||||
|
||||
def process_file(filename, treename, particle=None, max_events=None):
|
||||
tree = uproot.open(filename)[f"{treename}"]
|
||||
|
||||
|
|
@ -410,7 +401,6 @@ def process_file(filename, treename, particle=None, max_events=None):
|
|||
print(f"File {filename} particle {particle}, tree {treename}")
|
||||
return prepare_tree_data(tree, treename, particle, max_events=max_events)
|
||||
|
||||
|
||||
def power_fit_and_plot(x, y, label, color=None):
|
||||
|
||||
x = np.array(x)
|
||||
|
|
@ -488,11 +478,9 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
self.file = None
|
||||
print("\nATTENTION: Root file not found, continue without uproot functions or uproot file manually")
|
||||
|
||||
|
||||
#intro = "Interactive Shell Started. Type 'help' to see commands."
|
||||
prompt = ">> "
|
||||
|
||||
|
||||
def default(self, line):
|
||||
# Check if the command starts with our multi-word phrase
|
||||
if line.startswith("make table "):
|
||||
|
|
@ -761,7 +749,6 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
import readline
|
||||
import atexit
|
||||
import os
|
||||
import awkward as ak
|
||||
|
||||
histfile = os.path.expanduser("~/.uproot_shell_history")
|
||||
|
||||
|
|
@ -781,7 +768,6 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
"file": self.file,
|
||||
"uproot": uproot,
|
||||
"np": np,
|
||||
#"ak": ak,
|
||||
"plt": plt,
|
||||
"self": self
|
||||
}
|
||||
|
|
@ -820,7 +806,7 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
look up 'uproot' for more details"""
|
||||
self.run_command_line()
|
||||
|
||||
def do_energy_analysis(self, arg):
|
||||
def do_energy_analysis(self, arg): #geometry needs update
|
||||
|
||||
args = shlex.split(arg)
|
||||
|
||||
|
|
@ -841,7 +827,6 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
|
||||
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"
|
||||
|
|
@ -901,7 +886,7 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
global EA
|
||||
EA = energy_loss(particle, Ei, dA)
|
||||
global EC
|
||||
EC = energy_loss(particle, Ei, ClassSX3dC)
|
||||
EC = energy_loss(particle, Ei, dC)
|
||||
global Esx3
|
||||
Esx3 = energy_loss(particle, Ei, dsx3)
|
||||
global Eprop
|
||||
|
|
@ -910,25 +895,15 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
Elost = Ei - Esx3
|
||||
|
||||
print("Analysis complete")
|
||||
|
||||
print(f"Processed events: {len(Ei)}")
|
||||
|
||||
print(f"Anode average energy: {np.mean(EA):.3f} MeV")
|
||||
|
||||
print(f"Cathode average energy: {np.mean(EC):.3f} MeV")
|
||||
|
||||
print(f"sx3 average energy: {np.mean(Esx3):.3f} MeV")
|
||||
|
||||
print(f"Average total energy loClassSX3ss to sx3: {np.mean(Elost):.3f} MeV")
|
||||
|
||||
print(f"Maximum total energy loss to sx3: {np.max(Elost):.3f} MeV")
|
||||
|
||||
print(f"Minimum total energy loss to sx3: {np.min(Elost):.3f} MeV")
|
||||
|
||||
print(f"Proportion counter average energy difference: {np.mean(Eprop):.3f} MeV")
|
||||
|
||||
print(f"Maximum proportion counter energy difference: {np.max(Eprop):.3f} MeV")
|
||||
|
||||
print(f"Minimum proportion counter energy difference: {np.min(Eprop):.3f} MeV")
|
||||
|
||||
|
||||
|
|
@ -937,46 +912,30 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
print(f"Writing new tree to {output_filename}")
|
||||
|
||||
# Load ALL original branches
|
||||
all_data = self.tree.arrays(
|
||||
library="np",
|
||||
entry_stop=max_events
|
||||
)
|
||||
all_data = self.tree.arrays(library="np",entry_stop=max_events)
|
||||
|
||||
# Create full-length arrays initialized to NaN
|
||||
n_total = len(data["Tb"])
|
||||
|
||||
EA_full = np.full(n_total, np.nan)
|
||||
|
||||
EC_full = np.full(n_total, np.nan)
|
||||
|
||||
Esx3_full = np.full(n_total, np.nan)
|
||||
|
||||
Eprop_full = np.full(n_total, np.nan)
|
||||
|
||||
Elost_full = np.full(n_total, np.nan)
|
||||
|
||||
# Put values back into valid entries
|
||||
EA_full[mask] = EA
|
||||
|
||||
EC_full[mask] = EC
|
||||
|
||||
Esx3_full[mask] = Esx3
|
||||
|
||||
Eprop_full[mask] = Eprop
|
||||
|
||||
Elost_full[mask] = Elost
|
||||
|
||||
# Add new branches
|
||||
all_data["EA"] = EA_full
|
||||
|
||||
all_data["EC"] = EC_full
|
||||
|
||||
all_data["Esx3"] = Esx3_full
|
||||
|
||||
all_data["Eprop"] = Eprop_full
|
||||
|
||||
all_data["Elost"] = Elost_full
|
||||
|
||||
# Write new ROOT file as a classic TTree
|
||||
with uproot.recreate(output_filename) as fout:
|
||||
|
||||
|
|
@ -1010,17 +969,15 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
|
||||
data = prepare_tree_data(self.tree, treename, particle, max_events=max_events)
|
||||
|
||||
mask = data["thetab"] >= 0
|
||||
|
||||
Ei = data["Ei"][mask]
|
||||
sx3Z = data["sx3Z"][mask]
|
||||
EA = data["EA"][mask]
|
||||
EC = data["EC"][mask]
|
||||
Esx3 = data["Esx3"][mask]
|
||||
Elost = data["Elost"][mask]
|
||||
Eprop = data["Eprop"][mask]
|
||||
dA = data["dA"][mask]
|
||||
dC = data["dC"][mask]
|
||||
Ei = data["Ei"]
|
||||
sx3Z = data["sx3Z"]
|
||||
EA = data["EA"]
|
||||
EC = data["EC"]
|
||||
Esx3 = data["Esx3"]
|
||||
Elost = data["Elost"]
|
||||
Eprop = data["Eprop"]
|
||||
dA = data["dA"]
|
||||
dC = data["dC"]
|
||||
|
||||
update_plot_data(f"{particle}_{treename}_Ei", Ei)
|
||||
update_plot_data(f"{particle}_{treename}_sx3Z", sx3Z)
|
||||
|
|
@ -1170,7 +1127,7 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
|
||||
# Keep only events with thetab >= 45
|
||||
#thetab_mask = all_branches["thetab"] >= 0
|
||||
thetab_mask = mask
|
||||
|
||||
|
||||
for branch in branch_names:
|
||||
values = all_branches[branch]
|
||||
|
|
@ -1186,7 +1143,7 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
continue
|
||||
|
||||
# Apply the thetab cut
|
||||
values = values[thetab_mask]
|
||||
values = values
|
||||
|
||||
values = values[~np.isnan(values)]
|
||||
|
||||
|
|
@ -1212,10 +1169,6 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
|
||||
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)
|
||||
|
|
@ -1287,7 +1240,20 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
#tree2_name = args[3] if len(args) > 3 else 'tree1'
|
||||
outdir = "dual_plots"
|
||||
|
||||
data1 = process_file(os.path.join("..", "Armory", file1), "tree1")
|
||||
# load both trees for file1 and combine their arrays into a single dataset
|
||||
data1_tree1 = process_file(os.path.join("..", "Armory", file1), "tree1")
|
||||
data1_tree2 = process_file(os.path.join("..", "Armory", file1), "tree2")
|
||||
# concatenate matching array entries
|
||||
data1 = {"particle": f"{data1_tree1['particle']}_combined"}
|
||||
for key in data1_tree1:
|
||||
if key == "particle":
|
||||
continue
|
||||
try:
|
||||
data1[key] = np.concatenate([data1_tree1[key], data1_tree2[key]])
|
||||
except Exception:
|
||||
# fallback: prefer tree1 value if concat fails
|
||||
data1[key] = data1_tree1[key]
|
||||
|
||||
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'])}")
|
||||
|
|
@ -1295,27 +1261,10 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
|
||||
print(f"Saving plots to: {outdir}")
|
||||
|
||||
|
||||
#Overlay histogram: Elost
|
||||
plt.figure(figsize=(8,6))
|
||||
|
||||
plt.hist(
|
||||
data1["Elost"],
|
||||
bins=200,
|
||||
histtype='step',
|
||||
linewidth=2,
|
||||
density=True,
|
||||
label=data1["particle"]
|
||||
)
|
||||
|
||||
plt.hist(
|
||||
data2["Elost"],
|
||||
bins=200,
|
||||
histtype='step',
|
||||
linewidth=2,
|
||||
density=True,
|
||||
label=data2["particle"]
|
||||
)
|
||||
plt.hist(data1["Elost"],bins=200,histtype='step',linewidth=2,density=True,label=data1["particle"])
|
||||
plt.hist(data2["Elost"],bins=200,histtype='step',linewidth=2,density=True,label=data2["particle"])
|
||||
|
||||
plt.xlabel("Energy Loss (MeV)")
|
||||
plt.ylabel("Normalized Counts")
|
||||
|
|
@ -1329,53 +1278,27 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
|
||||
#Overlay histogram: sx3Z
|
||||
plt.figure(figsize=(8,6))
|
||||
|
||||
plt.hist(
|
||||
data1["sx3Z"],
|
||||
bins=150,
|
||||
histtype='step',
|
||||
linewidth=2,
|
||||
density=True,
|
||||
label=data1["particle"]
|
||||
)
|
||||
|
||||
plt.hist(
|
||||
data2["sx3Z"],
|
||||
bins=150,
|
||||
histtype='step',
|
||||
linewidth=2,
|
||||
density=True,
|
||||
label=data2["particle"]
|
||||
)
|
||||
|
||||
plt.hist(data1["sx3Z"],bins=150,histtype='step',linewidth=2,density=True,label=data1["particle"])
|
||||
plt.hist(data2["sx3Z"],bins=150,histtype='step',linewidth=2,density=True,label=data2["particle"])
|
||||
plt.xlabel("SX3 Z")
|
||||
plt.ylabel("Normalized Counts")
|
||||
plt.title("SX3 Position Comparison")
|
||||
plt.legend()
|
||||
plt.grid(True)
|
||||
plt.tight_layout()
|
||||
|
||||
plt.savefig(f"{outdir}/sx3Z_overlay.png", dpi=300)
|
||||
plt.show()
|
||||
|
||||
try:
|
||||
#Side-by-side 2D plots
|
||||
fig, axes = plt.subplots(1, 2, figsize=(14,6))
|
||||
|
||||
h1 = axes[0].hist2d(
|
||||
data1["sx3Z"],
|
||||
data1["Elost"],
|
||||
bins=200
|
||||
)
|
||||
h1 = axes[0].hist2d(data1["sx3Z"],data1["Elost"],bins=200)
|
||||
|
||||
axes[0].set_title(f'{data1["particle"]} Elost vs SX3')
|
||||
axes[0].set_xlabel("SX3 Z")
|
||||
axes[0].set_ylabel("Energy Loss (MeV)")
|
||||
|
||||
h2 = axes[1].hist2d(
|
||||
data2["sx3Z"],
|
||||
data2["Elost"],
|
||||
bins=200
|
||||
)
|
||||
h2 = axes[1].hist2d(data2["sx3Z"],data2["Elost"],bins=200)
|
||||
|
||||
axes[1].set_title(f'{data2["particle"]} Elost vs SX3')
|
||||
axes[1].set_xlabel("SX3 Z")
|
||||
|
|
@ -1388,47 +1311,48 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
|
||||
plt.savefig(f"{outdir}/Elost_vs_sx3_comparison.png", dpi=300)
|
||||
plt.show()
|
||||
except:
|
||||
print("Error with side-by-side plots")
|
||||
|
||||
#EA vs Esx3 overlay scatter
|
||||
plt.figure(figsize=(8,6))
|
||||
|
||||
plt.scatter(
|
||||
data1["EA"],
|
||||
data1["Esx3"],
|
||||
s=1,
|
||||
alpha=0.3,
|
||||
label=data1["particle"]
|
||||
)
|
||||
|
||||
plt.scatter(
|
||||
data2["EA"],
|
||||
data2["Esx3"],
|
||||
s=1,
|
||||
alpha=0.3,
|
||||
label=data2["particle"]
|
||||
)
|
||||
|
||||
plt.scatter(data1["EA"],data1["Esx3"],s=1,alpha=0.3,label=data1["particle"])
|
||||
plt.scatter(data2["EA"],data2["Esx3"],s=1,alpha=0.3,label=data2["particle"])
|
||||
plt.xlabel("EA (MeV)")
|
||||
plt.ylabel("Esx3 (MeV)")
|
||||
plt.title("Anode vs SX3 Energy")
|
||||
plt.legend()
|
||||
plt.grid(True)
|
||||
plt.tight_layout()
|
||||
|
||||
plt.savefig(f"{outdir}/EA_vs_Esx3_overlay.png", dpi=300)
|
||||
plt.show()
|
||||
|
||||
#PCE vs SiE
|
||||
|
||||
mask1 = data1["Esx3"] > 0
|
||||
mask2 = data2["Esx3"] > 0
|
||||
mask1 = data1["Esx3"] > 1
|
||||
mask2 = data2["Esx3"] > 1
|
||||
|
||||
thetab1 = np.deg2rad(data1["thetab"][mask1])
|
||||
thetab2 = np.deg2rad(data2["thetab"][mask2])
|
||||
|
||||
#theta smear (spatial uncertainty)
|
||||
if True:
|
||||
sigma = 0.2 * np.sqrt(np.maximum(thetab1, 0))
|
||||
|
||||
thetab1 += np.random.normal(0,sigma,len(thetab1))
|
||||
thetab1 += np.random.normal(0,0.02*np.sqrt(thetab1),len(thetab1))
|
||||
|
||||
sigma = 0.2 * np.sqrt(np.maximum(thetab2, 0))
|
||||
thetab2 += np.random.normal(0,sigma,len(thetab2))
|
||||
thetab2 += np.random.normal(0,0.02*np.sqrt(thetab2),len(thetab2))
|
||||
|
||||
combined_Esx3 = np.concatenate([
|
||||
(data1["Esx3"][mask1]),
|
||||
data2["Esx3"][mask2]
|
||||
])
|
||||
combined_Eprop = np.concatenate([
|
||||
data1["Eprop"][mask1] * 3,
|
||||
data2["Eprop"][mask2] * 2.7
|
||||
data1["Eprop"][mask1] * np.sin(thetab1) * 3,
|
||||
data2["Eprop"][mask2] * np.sin(thetab2) * 3
|
||||
])
|
||||
|
||||
combined_Esx3 = (
|
||||
|
|
@ -1436,7 +1360,8 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
+ np.random.normal(0,0.08,len(combined_Esx3))
|
||||
)
|
||||
|
||||
if True:
|
||||
#Artifical smear
|
||||
if False:
|
||||
sigma = 0.02 * np.sqrt(np.maximum(combined_Eprop, 0))
|
||||
|
||||
combined_Eprop += np.random.normal(
|
||||
|
|
@ -1451,13 +1376,9 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
len(combined_Eprop)
|
||||
)
|
||||
|
||||
mask = (
|
||||
np.isfinite(combined_Esx3)
|
||||
& np.isfinite(combined_Eprop)
|
||||
)
|
||||
|
||||
mask = (np.isfinite(combined_Esx3)& np.isfinite(combined_Eprop))
|
||||
combined_Esx3 = combined_Esx3[mask]
|
||||
combined_Eprop = combined_Eprop[mask]# * 3
|
||||
combined_Eprop = combined_Eprop[mask]
|
||||
|
||||
#combined_Eprop = combined_Eprop * .686
|
||||
|
||||
|
|
@ -1468,27 +1389,27 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
bins=200
|
||||
)
|
||||
plt.xlabel("SX3 Energy (MeV)")
|
||||
plt.ylabel("PCEnergy")
|
||||
plt.xlim(0,30)
|
||||
plt.ylim(0,0.45)
|
||||
plt.ylabel("PCEnergy x Sin(theta)")
|
||||
#plt.xlim(0,30)
|
||||
#plt.ylim(0,0.45)
|
||||
plt.colorbar(label="Counts")
|
||||
plt.tight_layout()
|
||||
plt.savefig(f"{outdir}/Eprop_vs_Esx3.png", dpi=300)
|
||||
plt.show()
|
||||
|
||||
plt.figure(figsize=(11,6), facecolor='white')
|
||||
plt.figure(figsize=(14,6), facecolor='white')
|
||||
plt.hist2d(
|
||||
combined_Esx3,
|
||||
combined_Eprop,
|
||||
bins=[500,500],
|
||||
#range=[[0,35],[0,0.6]],
|
||||
cmap=yellow_jet,
|
||||
cmin=1
|
||||
)
|
||||
cmin=1)
|
||||
plt.margins(0)
|
||||
plt.xlabel("SX3 Energy (MeV)")
|
||||
plt.ylabel("PCEnergy")
|
||||
plt.xlim(0,30)
|
||||
plt.ylim(0,0.45)
|
||||
plt.ylabel("PCEnergy x Sin(theta)")
|
||||
plt.xlim(0,35)
|
||||
plt.ylim(0,0.6)
|
||||
cbar = plt.colorbar()
|
||||
cbar.set_label("Counts")
|
||||
plt.minorticks_on()
|
||||
|
|
@ -1513,14 +1434,14 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
mask = data1["Esx3"] > 0
|
||||
power_fit_and_plot(
|
||||
data1["Esx3"][mask],
|
||||
data1["Eprop"][mask] * .686,
|
||||
data1["Eprop"][mask],
|
||||
label=data1["particle"],
|
||||
color='red'
|
||||
)
|
||||
mask2 = data2["Esx3"] > 0
|
||||
power_fit_and_plot(
|
||||
data2["Esx3"][mask2],
|
||||
data2["Eprop"][mask2] * .686,
|
||||
data2["Eprop"][mask2],
|
||||
label=data2["particle"],
|
||||
color='cyan'
|
||||
)
|
||||
|
|
@ -1539,11 +1460,8 @@ class MyInteractiveApp(cmd.Cmd):
|
|||
plt.savefig(f"{outdir}/Combined_Eprop_vs_Esx3_fit.png",dpi=300)
|
||||
plt.show()
|
||||
|
||||
|
||||
|
||||
print("Dual plotting complete.")
|
||||
|
||||
|
||||
#exec(open("PCEnergyAnalysis.py").read())
|
||||
|
||||
if __name__ == "__main__":
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user