rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

anasen cpp updates

Browse Source
This commit is contained in:
James Szalkie 2026-08-13 12:42:50 -04:00
parent b5add710e5
commit cb22da5f57
2 changed files with 224 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

88
Armory/anasenMS.cpp
View File

@ -15,6 +15,7 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <set> #include <set>
#include <tuple>
#include "TLegend.h" #include "TLegend.h"
#include "TH1D.h" #include "TH1D.h"
#include "TObjArray.h" #include "TObjArray.h"
@ -42,10 +43,38 @@ bool IsDeadCathode(int id){
} }
bool IsDeadSX3(int id){ bool IsDeadSX3(int id){
static std::set<int> dead = {}; // add dead SX3 IDs here, 0-23 1,7,9,3 static std::set<int> dead = {}; //{0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23}; // add dead SX3 IDs here, 0-23 1,7,9,3
return dead.count(id); return dead.count(id);
} }
bool IsDeadSX3FrontUpChannel(int sx3ID, int chUp){
static std::set<std::pair<int, int>> dead = {
// {sx3ID, front-up-channel}
};
return dead.count({sx3ID, chUp});
}
bool IsDeadSX3FrontDnChannel(int sx3ID, int chDn){
static std::set<std::pair<int, int>> dead = {
// {sx3ID, front-down-channel}
};
return dead.count({sx3ID, chDn});
}
bool IsDeadSX3BackChannel(int sx3ID, int chBk){
static std::set<std::pair<int, int>> dead = {
//{9, 10}
// {sx3ID, back-channel}
};
return dead.count({sx3ID, chBk});
}
bool IsDeadSX3ChannelCombo(int sx3ID, int chUp, int chDn, int chBk){
return IsDeadSX3FrontUpChannel(sx3ID, chUp)
|| IsDeadSX3FrontDnChannel(sx3ID, chDn)
|| IsDeadSX3BackChannel(sx3ID, chBk);
}
// Simulate sequential two-body decay of an unstable parent in its rest frame. // Simulate sequential two-body decay of an unstable parent in its rest frame.
TLorentzVector SimulateSequentialDecay(const TLorentzVector &parent, TLorentzVector SimulateSequentialDecay(const TLorentzVector &parent,
int daughterA, int daughterZ, int daughterA, int daughterZ,
@ -98,29 +127,33 @@ int main(int argc, char **argv){
TransferReaction transfer; TransferReaction transfer;
//To set beam energy loss, use energy loss app, and create table with target isotope, set Initial beam energy as max energy //To set beam energy loss, use energy loss app, and create table with target isotope, set Initial beam energy as max energy
transfer.SetA(18, 9, 0); // 18Ne projectile transfer.SetA(18, 10 0); // 22Mg projectile
//TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_Na-21.dat"); TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_Ne-18.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(21, 11); // 21Na* heavy product transfer.SetB(21, 11); // 30Si* heavy product
const ReactionConfig reactionConfig = transfer.GetRectionConfig(); const ReactionConfig reactionConfig = transfer.GetRectionConfig();
const double beamA = reactionConfig.beamA; // mass number of 14N beam const double beamA = reactionConfig.beamA; // mass number of 14N beam
const double beamE = 3; // maximum beam energy in MeV/u //const double beamE = 72 / beamA; // beam energy in MeV
// Excited state lists (projectile and heavy-product excitation states)
std::vector<float> ExAList = {0}; // Beam excited energy
std::vector<float> ExList = {0, 0.3, 1.7, 2.4, 2.8}; // Heavy product excited energy
const int kMBeam = reactionConfig.beamA; // mass number of beam const int kMBeam = reactionConfig.beamA; // mass number of beam
const int kMTarget = reactionConfig.targetA; // mass number of target const int kMTarget = reactionConfig.targetA; // mass number of target
const int kMLight = reactionConfig.recoilLightA; // mass number of light ejectile const int kMLight = reactionConfig.recoilLightA; // mass number of light ejectile
const int kMHeavy = reactionConfig.recoilHeavyA; // mass number of heavy product const int kMHeavy = reactionConfig.recoilHeavyA; // mass number of heavy product
const int kZBeam = reactionConfig.beamZ; // atomic number of beam
const int kZTarget = reactionConfig.targetZ; // atomic number of target
const int kZLight = reactionConfig.recoilLightZ; // atomic number of light ejectile
const int kZHeavy = reactionConfig.recoilHeavyZ; // atomic number of heavy product
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;
const int decayDaughterZ = 10; const int decayDaughterZ = 10;
const int decayEjectA = 1; const int decayEjectA = 1;
const int decayEjectZ = 1; const int decayEjectZ = 1;
// Excited state lists (projectile and heavy-product excitation states)
std::vector<float> ExAList = {0}; // Beam excited energy
std::vector<float> ExList = {0, .3, 1.7, 2.4, 2.8, 3.5, 3.9, 4, 4.3, 4.5}; // Heavy product excited energy
// 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
@ -188,6 +221,10 @@ int main(int argc, char **argv){
int MTargetOut; int MTargetOut;
int MLightOut; int MLightOut;
int MHeavyOut; int MHeavyOut;
int ZBeamOut;
int ZTargetOut;
int ZLightOut;
int ZHeavyOut;
tree1->Branch("beamKEA", &KEA, "beamKEA/D"); tree1->Branch("beamKEA", &KEA, "beamKEA/D");
tree2->Branch("beamKEA", &KEA2, "beamKEA/D"); tree2->Branch("beamKEA", &KEA2, "beamKEA/D");
tree1->Branch("beamPath_cm", &beamPath_cm, "beamPath_cm/D"); tree1->Branch("beamPath_cm", &beamPath_cm, "beamPath_cm/D");
@ -200,6 +237,10 @@ int main(int argc, char **argv){
tree1->Branch("MTarget", &MTargetOut, "MTarget/I"); tree1->Branch("MTarget", &MTargetOut, "MTarget/I");
tree1->Branch("MLight", &MLightOut, "MLight/I"); tree1->Branch("MLight", &MLightOut, "MLight/I");
tree1->Branch("MHeavy", &MHeavyOut, "MHeavy/I"); tree1->Branch("MHeavy", &MHeavyOut, "MHeavy/I");
tree1->Branch("ZBeam", &ZBeamOut, "ZBeam/I");
tree1->Branch("ZTarget", &ZTargetOut, "ZTarget/I");
tree1->Branch("ZLight", &ZLightOut, "ZLight/I");
tree1->Branch("ZHeavy", &ZHeavyOut, "ZHeavy/I");
tree2->Branch("MBeam", &MBeamOut, "MBeam/I"); tree2->Branch("MBeam", &MBeamOut, "MBeam/I");
tree2->Branch("MTarget", &MTargetOut, "MTarget/I"); tree2->Branch("MTarget", &MTargetOut, "MTarget/I");
tree2->Branch("MLight", &MLightOut, "MLight/I"); tree2->Branch("MLight", &MLightOut, "MLight/I");
@ -388,12 +429,10 @@ int main(int argc, char **argv){
// compute beam energy at the event vertex from the gas path length // compute beam energy at the event vertex from the gas path length
beamPath_cm = TVector3(vertexZ - beamEntranceZ, vertexX, vertexY).Mag() * 0.1; beamPath_cm = TVector3(vertexZ - beamEntranceZ, vertexX, vertexY).Mag() * 0.1;
if( beamPath_cm < 0 ) beamPath_cm = 0; if( beamPath_cm < 0 ) beamPath_cm = 0;
//beamEnergy = elossBeam->Eval(beamPath_cm); // MeV beamEnergy = elossBeam->Eval(beamPath_cm); // MeV
//beamEnergyLoss = elossBeam->Eval(0.0) - beamEnergy; beamEnergyLoss = elossBeam->Eval(0.0) - beamEnergy;
//KEA = beamEnergy / beamA; KEA = beamEnergy / beamA;
KEA = gRandom->Uniform(0, beamE); //KEA = gRandom->Uniform(0, beamE);
beamEnergy = KEA * beamA;
beamEnergyLoss = 0;
transfer.SetIncidentEnergyAngle(KEA, 0, 0); transfer.SetIncidentEnergyAngle(KEA, 0, 0);
transfer.CalReactionConstant(); transfer.CalReactionConstant();
@ -479,6 +518,7 @@ int main(int argc, char **argv){
sx3Up = sx3->GetChUp(); sx3Up = sx3->GetChUp();
sx3Dn = sx3->GetChDn(); sx3Dn = sx3->GetChDn();
sx3Bk = sx3->GetChBk(); sx3Bk = sx3->GetChBk();
if(IsDeadSX3ChannelCombo(sx3ID, sx3Up, sx3Dn, sx3Bk)) continue;
sx3ZFrac = sx3->GetZFrac(); sx3ZFrac = sx3->GetZFrac();
// apply intrinsic detector resolution to true SX3 hit position // apply intrinsic detector resolution to true SX3 hit position
@ -555,6 +595,23 @@ int main(int argc, char **argv){
sx3Up2 = sx3->GetChUp(); sx3Up2 = sx3->GetChUp();
sx3Dn2 = sx3->GetChDn(); sx3Dn2 = sx3->GetChDn();
sx3Bk2 = sx3->GetChBk(); sx3Bk2 = sx3->GetChBk();
if(IsDeadSX3ChannelCombo(sx3ID2, sx3Up2, sx3Dn2, sx3Bk2)){
sx3ID2 = -1;
sx3Up2 = -1;
sx3Dn2 = -1;
sx3Bk2 = -1;
sx3ZFrac2 = TMath::QuietNaN();
sx3X2 = TMath::QuietNaN();
sx3Y2 = TMath::QuietNaN();
sx3Z2 = TMath::QuietNaN();
anodeDist2[0] = TMath::QuietNaN();
cathodeDist2[0] = TMath::QuietNaN();
reTheta2 = TMath::QuietNaN();
rePhi2 = TMath::QuietNaN();
reTheta12 = TMath::QuietNaN();
rePhi12 = TMath::QuietNaN();
z02 = TMath::QuietNaN();
} else {
sx3ZFrac2 = sx3->GetZFrac(); sx3ZFrac2 = sx3->GetZFrac();
TVector3 hitPos2 = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L); TVector3 hitPos2 = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L);
sx3X2 = hitPos2.X(); sx3X2 = hitPos2.X();
@ -567,6 +624,7 @@ int main(int argc, char **argv){
reTheta12 = pw->GetTrackTheta() * TMath::RadToDeg(); reTheta12 = pw->GetTrackTheta() * TMath::RadToDeg();
rePhi12 = pw->GetTrackPhi() * TMath::RadToDeg(); rePhi12 = pw->GetTrackPhi() * TMath::RadToDeg();
z02 = pw->GetZ0(); z02 = pw->GetZ0();
}
} }
KEA2 = KEA; KEA2 = KEA;

176
ELoss/PCEnergyAnalysis.py
View File

@ -78,6 +78,45 @@ SX3_SI_THICKNESS_CM = 0.1
QQQ_SI_THICKNESS_CM = 0.1 QQQ_SI_THICKNESS_CM = 0.1
SX3_THETA_MIN_DEG = 0.0 SX3_THETA_MIN_DEG = 0.0
# Optional dead-channel masks for SX3.
# Format for channel masks is: (sx3ID, channel)
DEAD_SX3_IDS = set() #Ex {9}
DEAD_SX3_FRONT_UP = set() # Ex {(9, 3)}
DEAD_SX3_FRONT_DN = set() # Ex {(9, 3)}
DEAD_SX3_BACK = set() # Ex {(9, 3)}
def _normalize_channel_array(values):
arr = np.asarray(values)
if arr.ndim > 1:
return arr[:, 0]
return arr
def build_sx3_live_mask(sx3_id, sx3_up=None, sx3_dn=None, sx3_bk=None):
sx3_id = _normalize_channel_array(sx3_id)
live_mask = np.ones(sx3_id.shape[0], dtype=bool)
if DEAD_SX3_IDS:
live_mask &= ~np.isin(sx3_id, list(DEAD_SX3_IDS))
if sx3_up is not None and DEAD_SX3_FRONT_UP:
sx3_up = _normalize_channel_array(sx3_up)
for det_id, ch in DEAD_SX3_FRONT_UP:
live_mask &= ~((sx3_id == det_id) & (sx3_up == ch))
if sx3_dn is not None and DEAD_SX3_FRONT_DN:
sx3_dn = _normalize_channel_array(sx3_dn)
for det_id, ch in DEAD_SX3_FRONT_DN:
live_mask &= ~((sx3_id == det_id) & (sx3_dn == ch))
if sx3_bk is not None and DEAD_SX3_BACK:
sx3_bk = _normalize_channel_array(sx3_bk)
for det_id, ch in DEAD_SX3_BACK:
live_mask &= ~((sx3_id == det_id) & (sx3_bk == ch))
return live_mask
def get_loss_table_path(medium, particle_label): def get_loss_table_path(medium, particle_label):
script_dir = os.path.dirname(os.path.abspath(__file__)) script_dir = os.path.dirname(os.path.abspath(__file__))
@ -451,15 +490,24 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
vertex_x = data["vX"] vertex_x = data["vX"]
vertex_y = data["vY"] vertex_y = data["vY"]
vertex_z = data["vZ"] vertex_z = data["vZ"]
beamEnergy = data["beamEnergy"]
sx3_theta_rad = np.radians(data["thetab"]) sx3_theta_rad = np.radians(data["thetab"])
sx3_phi_rad = np.radians(data["phib"]) sx3_phi_rad = np.radians(data["phib"])
sx3_hit_x = data["sx3X"] sx3_hit_x = data["sx3X"]
sx3_hit_y = data["sx3Y"] sx3_hit_y = data["sx3Y"]
sx3_hit_z = data["sx3Z"] sx3_hit_z = data["sx3Z"]
sx3_id = data["sx3ID"] if "sx3ID" in data else np.full_like(sx3_hit_x, -1)
sx3_up = data["sx3Up"] if "sx3Up" in data else None
sx3_dn = data["sx3Dn"] if "sx3Dn" in data else None
sx3_bk = data["sx3Bk"] if "sx3Bk" in data else None
sx3_angle_mask = data["thetab"] > SX3_THETA_MIN_DEG sx3_angle_mask = data["thetab"] > SX3_THETA_MIN_DEG
sx3_valid_mask = ~np.isnan(sx3_hit_x) & ~np.isnan(sx3_hit_y) & ~np.isnan(sx3_hit_z) sx3_valid_mask = ~np.isnan(sx3_hit_x) & ~np.isnan(sx3_hit_y) & ~np.isnan(sx3_hit_z)
sx3_event_mask = sx3_angle_mask & sx3_valid_mask sx3_live_mask = build_sx3_live_mask(sx3_id, sx3_up, sx3_dn, sx3_bk)
sx3_event_mask = sx3_angle_mask & sx3_valid_mask & sx3_live_mask
n_rejected_dead = int(np.sum(sx3_angle_mask & sx3_valid_mask & ~sx3_live_mask))
if n_rejected_dead > 0:
print(f"Rejected {n_rejected_dead} SX3 events due to dead-channel mask")
qqq_hit_x = data["qqqX"] qqq_hit_x = data["qqqX"]
qqq_hit_y = data["qqqY"] qqq_hit_y = data["qqqY"]
@ -632,7 +680,12 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
"Efinal": Efinal, "Efinal": Efinal,
"lsx3": sx3_silicon_path_cm, "lsx3": sx3_silicon_path_cm,
"lqqq": qqq_silicon_path_cm, "lqqq": qqq_silicon_path_cm,
"EdetQ": EdetQ "EdetQ": EdetQ,
"beamEnergy": beamEnergy,
"vZ": vertex_z,
"beamEnergyLoss": data["beamEnergyLoss"],
"Ex": data["Ex"],
"vZsx3": sx3_vertex_z
} }
def infer_particle_from_filename(filename): def infer_particle_from_filename(filename):
@ -1105,6 +1158,10 @@ class MyInteractiveApp(cmd.Cmd):
"vZ", "vZ",
"sx3X", "sx3X",
"sx3Y", "sx3Y",
"sx3ID",
"sx3Up",
"sx3Dn",
"sx3Bk",
"sx3XExit", "sx3XExit",
"sx3YExit", "sx3YExit",
"sx3ZExit", "sx3ZExit",
@ -1140,7 +1197,13 @@ class MyInteractiveApp(cmd.Cmd):
if treename == 'tree1': if treename == 'tree1':
sx3_hit_x = data["sx3X"] sx3_hit_x = data["sx3X"]
sx3_hit_y = data["sx3Y"] sx3_hit_y = data["sx3Y"]
event_mask = ~np.isnan(sx3_hit_x) & ~np.isnan(sx3_hit_y) & ~np.isnan(sx3Z) sx3_live_mask = build_sx3_live_mask(
data["sx3ID"],
data["sx3Up"],
data["sx3Dn"],
data["sx3Bk"],
)
event_mask = ~np.isnan(sx3_hit_x) & ~np.isnan(sx3_hit_y) & ~np.isnan(sx3Z) & sx3_live_mask
else: else:
event_mask = ~np.isnan(initial_energy) & ~np.isnan(polar_angle_rad) event_mask = ~np.isnan(initial_energy) & ~np.isnan(polar_angle_rad)
@ -1297,6 +1360,11 @@ class MyInteractiveApp(cmd.Cmd):
Elostqqq = data["Elostqqq"] Elostqqq = data["Elostqqq"]
Edet = data["Edet"] Edet = data["Edet"]
EdetQ = data["EdetQ"] EdetQ = data["EdetQ"]
beamEnergy = data["beamEnergy"]
vZ = data["vZ"]
beamEnergyLoss = data["beamEnergyLoss"]
Ex = data["Ex"]
vZsx3 = data["vZsx3"]
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)
@ -1324,7 +1392,7 @@ class MyInteractiveApp(cmd.Cmd):
sx3_theta_plot_mask = np.isfinite(thetab) & (thetab > SX3_THETA_MIN_DEG) sx3_theta_plot_mask = np.isfinite(thetab) & (thetab > SX3_THETA_MIN_DEG)
qqq_theta_plot_mask = np.isfinite(thetabqqq) & (thetabqqq > SX3_THETA_MIN_DEG) qqq_theta_plot_mask = np.isfinite(thetabqqq) & (thetabqqq > SX3_THETA_MIN_DEG)
if x.size > 0 and False: if x.size > 0 and True:
fig = plt.figure(figsize=(8,6)) fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection='3d') ax = fig.add_subplot(111, projection='3d')
@ -1422,8 +1490,6 @@ class MyInteractiveApp(cmd.Cmd):
plt.xlabel("Eqqq (MeV)") plt.xlabel("Eqqq (MeV)")
plt.title(f"{particle} ({treename}) Energy QQQ vs Elost * Theta") plt.title(f"{particle} ({treename}) Energy QQQ vs Elost * Theta")
plt.colorbar(label="Counts") plt.colorbar(label="Counts")
#plt.xlim(0,30)
#plt.ylim(0,0.45)
plt.tight_layout() plt.tight_layout()
plt.savefig(f"{base}/Eqqq_vs_Elostxsintheta.png", dpi=300) plt.savefig(f"{base}/Eqqq_vs_Elostxsintheta.png", dpi=300)
plt.show() plt.show()
@ -1436,8 +1502,6 @@ class MyInteractiveApp(cmd.Cmd):
plt.ylabel("Esx3") plt.ylabel("Esx3")
plt.title(f"{particle} ({treename}) sx3 Energy vs Theta") plt.title(f"{particle} ({treename}) sx3 Energy vs Theta")
plt.colorbar(label="Counts") plt.colorbar(label="Counts")
#plt.xlim(0,30)
#plt.ylim(0,0.45)
plt.tight_layout() plt.tight_layout()
plt.savefig(f"{base}/sx3E_vs_theta.png", dpi=300) plt.savefig(f"{base}/sx3E_vs_theta.png", dpi=300)
plt.show() plt.show()
@ -1650,8 +1714,6 @@ class MyInteractiveApp(cmd.Cmd):
plt.xlabel("SX3 Energy Detected (MeV)") plt.xlabel("SX3 Energy Detected (MeV)")
plt.title(f"{particle} ({treename}) Energy Propagation Difference vs SX3 Energy Detected") plt.title(f"{particle} ({treename}) Energy Propagation Difference vs SX3 Energy Detected")
plt.colorbar(label="Counts") plt.colorbar(label="Counts")
#plt.xlim(0,30)
#plt.ylim(0,.5)
plt.tight_layout() plt.tight_layout()
plt.savefig(f"{base}/Eprop_vs_Edet_log.png", dpi=300) plt.savefig(f"{base}/Eprop_vs_Edet_log.png", dpi=300)
plt.show() plt.show()
@ -1672,8 +1734,6 @@ class MyInteractiveApp(cmd.Cmd):
plt.xlabel("Eqqq detected (MeV)") plt.xlabel("Eqqq detected (MeV)")
plt.title(f"{particle} ({treename}) EDet QQQ vs PCEnergy * sin(theta)") plt.title(f"{particle} ({treename}) EDet QQQ vs PCEnergy * sin(theta)")
plt.colorbar(label="Counts") plt.colorbar(label="Counts")
#plt.xlim(0,30)
#plt.ylim(0,0.45)
plt.tight_layout() plt.tight_layout()
plt.savefig(f"{base}/Eprop_vs_EdetQ.png", dpi=300) plt.savefig(f"{base}/Eprop_vs_EdetQ.png", dpi=300)
plt.show() plt.show()
@ -1691,12 +1751,76 @@ class MyInteractiveApp(cmd.Cmd):
plt.xlabel("Eqqq detected (MeV)") plt.xlabel("Eqqq detected (MeV)")
plt.title(f"{particle} ({treename}) EDet QQQ vs PCEnergy * sin(theta)") plt.title(f"{particle} ({treename}) EDet QQQ vs PCEnergy * sin(theta)")
plt.colorbar(label="Counts") plt.colorbar(label="Counts")
#plt.xlim(0,30)
#plt.ylim(0,0.45)
plt.tight_layout() plt.tight_layout()
plt.savefig(f"{base}/Eprop_vs_EdetQ_log.png", dpi=300) plt.savefig(f"{base}/Eprop_vs_EdetQ_log.png", dpi=300)
plt.show() plt.show()
mask1 = (beamEnergy > 0) & ~np.isnan(vZ) & ~np.isnan(beamEnergy)
plt.figure(figsize=(7,6))
plt.hist2d(vZ[mask1], beamEnergy[mask1], bins=200)
plt.ylabel("Beam Energy")
plt.xlabel("Z")
plt.title(f"{particle} ({treename}) Beam Energy vs. Z")
plt.colorbar(label="Counts")
plt.tight_layout()
plt.savefig(f"{base}/EBeam_vs_Z", dpi=300)
plt.show()
mask1 = (beamEnergy > 0) & ~np.isnan(vZ) & ~np.isnan(beamEnergy)
plt.figure(figsize=(7,6))
plt.hist2d(vZ[mask1], beamEnergy[mask1], bins=200)
plt.ylabel("Beam Energy")
plt.xlabel("Z")
plt.title(f"{particle} ({treename}) Beam Energy vs. Z")
plt.colorbar(label="Counts")
plt.tight_layout()
plt.savefig(f"{base}/EBeam_vs_Z", dpi=300)
plt.show()
mask1 = (beamEnergy > 0) & ~np.isnan(vZ) & ~np.isnan(beamEnergy)
plt.figure(figsize=(7,6))
plt.hist2d(vZ[mask1], beamEnergy[mask1], bins=200)
plt.ylabel("Beam Energy")
plt.xlabel("Z")
plt.title(f"{particle} ({treename}) Beam Energy vs. Z")
plt.colorbar(label="Counts")
plt.tight_layout()
plt.savefig(f"{base}/EBeam_vs_Z", dpi=300)
plt.show()
mask1 = (beamEnergyLoss > 0) & ~np.isnan(vZ) & ~np.isnan(beamEnergy)
plt.figure(figsize=(7,6))
plt.hist2d(vZ[mask1], beamEnergyLoss[mask1], bins=200)
plt.ylabel("Beam Energy Loss")
plt.xlabel("Z")
plt.title(f"{particle} ({treename}) Beam Energy vs. Z")
plt.colorbar(label="Counts")
plt.tight_layout()
plt.savefig(f"{base}/EBeamLoss_vs_Z", dpi=300)
plt.show()
plt.figure(figsize=(7,6))
plt.hist2d(vZ, Ex, bins=200)
plt.ylabel("Excitation Energy")
plt.xlabel("Z")
#plt.ylim(0, 5)
plt.title(f"{particle} ({treename}) Excitation Energy vs. Z")
plt.colorbar(label="Counts")
plt.tight_layout()
plt.savefig(f"{base}/Ex_vs_Z", dpi=300)
plt.show()
#mask1 = ~np.isnan(Elost) & ~np.isnan(vZ)
plt.figure(figsize=(7,6))
plt.hist2d(vZsx3, Elost, bins=200)
plt.ylabel("Energy loss sx3")
plt.xlabel("Z")
#plt.ylim(0, 5)
plt.title(f"{particle} ({treename}) Energy loss vs. Z")
plt.colorbar(label="Counts")
plt.tight_layout()
plt.savefig(f"{base}/Eloss_vs_Z", dpi=300)
plt.show()
branch_names = [] branch_names = []
for key in self.tree.keys(): for key in self.tree.keys():
@ -1733,21 +1857,23 @@ class MyInteractiveApp(cmd.Cmd):
if values.size == 0: if values.size == 0:
print(f"Skipping empty branch: {branch}") print(f"Skipping empty branch: {branch}")
continue continue
plt.figure(figsize=(7,5)) plt.figure(figsize=(7,5))
try: try:
plt.hist(values, bins=100) 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()
except: except:
plt.hist(values, bins=10) print(f"Can not print branch {branch}")
plt.xlabel(branch) continue
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: else:
print("No branches found to histogram.") print("No branches found to histogram.")