Skip to content

Commit

Permalink
Merge branch 'main' into mdolce_ke_3d_proj_evd
Browse files Browse the repository at this point in the history
  • Loading branch information
@mdolce8
mdolce8 authored Sep 19, 2024
2 parents 6457c64 + 1d98a94 commit 14f8a6b
Show file tree
Hide file tree
Showing 41 changed files with 2,830 additions and 797 deletions.
96 changes: 60 additions & 36 deletions app/ConvertToTMSTree.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include "TTree.h"
#include "TGeoManager.h"
#include "TStopwatch.h"
#include "TParameter.h"

// EDepSim includes
#include "EDepSim/TG4Event.h"
Expand Down Expand Up @@ -82,23 +83,44 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) {
// Do we overlay events
bool Overlay = false;
// How many events do we want to overlay?
int nOverlays = 3;
int nOverlays = 50;
// The vector carrying our events that we want to overlay
std::vector<TMS_Event> overlay_events;

bool NerscOverlay = false;
TParameter<double>* spillPeriod_s = (TParameter<double>*)input->Get("spillPeriod_s");
if (spillPeriod_s != NULL) NerscOverlay = true;
double SpillPeriod = 0;
if (NerscOverlay) {
std::cout<<"Combining spills"<<std::endl;
SpillPeriod = spillPeriod_s->GetVal() * 1e9; // convert to ns
std::cout<<"Found spillSeriod_s of "<<SpillPeriod<<"ns"<<std::endl;
if (SpillPeriod < 1e7 || SpillPeriod > 1e11) {
std::cout<<"Fatal: Found spillSeriod_s that is unusually high or low. Expecting something like ~1.2e9 ns"<<std::endl;
exit(1);
}
TMS_Manager::GetInstance().Set_Nersc_Spill_Period(SpillPeriod);
}
int current_spill_number = 0;

for (; i < N_entries; ++i) {
if (N_entries <= 10 || i % (N_entries/10) == 0) {
std::cout << "Processed " << i << "/" << N_entries << " (" << double(i)*100./N_entries << "%)" << std::endl;
}

events->GetEntry(i);
// todo, gRoo has a different indexing than events with overlay
if (gRoo)
gRoo->GetEntry(i);

if (N_entries <= 10 || i % (N_entries/10) == 0) {
std::cout << "Processed " << i << "/" << N_entries << " (" << double(i)*100./N_entries << "%)" << std::endl;
}

// Todo: This should no longer be needed when this bug is fixed in the spill builder
// https://github.com/DUNE/2x2_sim/issues/54
event->EventId = i;
//if (event->Primaries.size() > 0)
// std::cout<<"Entry "<<i<<", interaction number of vtx 0: "<<event->Primaries[0].GetInteractionNumber()<<", vs event.EventId "<<event->EventId<<std::endl;

// Make a TMS event
TMS_Event tms_event = TMS_Event(*event);

// Fill up truth information from the GRooTracker object
if (gRoo){
tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, EvtVtx);
Expand All @@ -109,13 +131,38 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) {
overlay_events.push_back(tms_event);
continue;
}

// Keep filling up the vector if within spill
if (NerscOverlay) {
double next_spill_time = (current_spill_number + 0.5) * TMS_Manager::GetInstance().Get_Nersc_Spill_Period();
double current_spill_time = event->Primaries.begin()->Position.T();
// Check that this neutrino is within spill, but not last event
if (current_spill_time < next_spill_time && i != N_entries - 1) {
overlay_events.push_back(tms_event);
continue;
}
}

// Add event information and truth from another event
if (Overlay && i % nOverlays == 0) {
if (overlay_events.size() > 0) {
// Now loop over previous events
for (auto &event : overlay_events) tms_event.AddEvent(event);
std::cout<<"Overlaying "<<overlay_events.size()<<" events"<<std::endl;
// The first event should be the starting point so reverse it
std::reverse(overlay_events.begin(), overlay_events.end());
TMS_Event last_event = overlay_events.back();
overlay_events.pop_back();
for (auto &event : overlay_events) last_event.AddEvent(event);
overlay_events.clear();
// Now add this event as the first event in the next set
overlay_events.push_back(tms_event);
if (NerscOverlay) current_spill_number += 1;
// ... and make this event the combined spill called "last_event"
tms_event = last_event;
}

// Apply the det sim now, after overlaying events
// This doesn't work right now
tms_event.ApplyReconstructionEffects();

// Dump information
//tms_event.Print();
Expand All @@ -127,6 +174,7 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) {
TMS_ReadoutTreeWriter::GetWriter().Fill(tms_event);

int nslices = TMS_TimeSlicer::GetSlicer().RunTimeSlicer(tms_event);
std::cout<<"Sliced event "<<i<<" into "<<nslices<<" slices"<<std::endl;

// Check if this is not pileup
if (gRoo && event->Primaries.size() == 1 && tms_event.GetNVertices() == 1) {
Expand Down Expand Up @@ -164,40 +212,16 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) {
double visible_energy_from_vertex = map[primary_vertex_id];
tms_event_slice.SetTotalVisibleEnergyFromVertex(visible_energy_from_vertex);

if ((unsigned long) primary_vertex_id >= event->Primaries.size())
std::cout<<"Warning: primary_vertex_id "<<primary_vertex_id<<" is above event->Primaries.size "<<event->Primaries.size()<<std::endl;
else {
// Now set the remaining information from the gRoo tchain.
auto primary_vertex = event->Primaries[primary_vertex_id];
int interaction_number = primary_vertex.GetInteractionNumber();
gRoo->GetEntry(interaction_number);
tms_event_slice.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, EvtVtx);
// And the lepton info
int lepton_index = -1;
int current_index = 0;
for (auto particle : primary_vertex.Particles) {
int pdg = std::abs(particle.GetPDGCode());
if (pdg >= 11 && pdg <= 16) {
lepton_index = current_index;
break;
}
current_index += 1;
}
if (lepton_index >= 0) {
auto lepton = primary_vertex.Particles[lepton_index];
int lepton_pdg = lepton.GetPDGCode();
auto lepton_position = primary_vertex.GetPosition();
auto lepton_momentum = lepton.GetMomentum();
tms_event_slice.FillTrueLeptonInfo(lepton_pdg, lepton_position, lepton_momentum);
}
}
gRoo->GetEntry(primary_vertex_id);
tms_event_slice.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, EvtVtx);
tms_event_slice.SetLeptonInfoUsingVertexID(primary_vertex_id);
}
}

event_counter += 1;
int spill_number = i;
tms_event_slice.SetSpillNumber(spill_number);

// Try finding some tracks
TMS_TrackFinder::GetFinder().FindTracks(tms_event_slice);

Expand Down
13 changes: 10 additions & 3 deletions config/TMS_Default_Config.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,10 @@
MinPoints = 2
# How far away can a bar/plane combination be to include in a cluster
Epsilon = 2.5
# Amount of neighbouring hits a hit needs to have for pre-clustering (DBSCAN)
PreDBNeighbours = 9
# Maximal distance to hits to count as neighbours for pre-clustering (DBSCAN)
PreDBDistance = 2.5

[Recon.Hough]
# The highest number of Hough transforms we can run after each other
Expand All @@ -42,7 +46,7 @@
# dist = sqrt(xplanes*xplanes + zplanes*zplanes)
MinDist = 4.0
# Do we run DBSCAN clustering before the Hough transform and then run Hough on each cluster?
FirstCluster = false
FirstCluster = true #false
# Do we merge adjacent tracks afterwards?
MergeTracks = true
# Do we run AStar afterwards on start and end points
Expand All @@ -63,11 +67,11 @@

[Recon.TrackMatch3D]
# How many plane layers can the ends/starts of the two simple tracks be away from each other?
PlaneLimit = 3
PlaneLimit = 2
# How many bars can the ends/starts of the two simple tracks be away from each other?
BarLimit = 12
# How many ns can the ends/starts of the two simple tracks be away from each other?
TimeLimit = 30
TimeLimit = 16 #30
# How many ns can the ends/starts of a X simple track be away from U simple track?
XTimeLimit = 15 # (~10 ns time difference for 3.1m long bars both, plus a bit extra)
# What is the 'anchor point' of the scintillator bars in y [mm]?
Expand All @@ -86,6 +90,9 @@
#CostMetric = "DetectorNotZ"
IsGreedy = false

[Recon.Kalman]
Run = true # Whether we run the Kalman filter or not

[Recon.Stopping]
nLastHits = 5 # The number of hits at the end of track to look at for stopping
EnergyCut = 2.5 # The cut in MeV to see if a track has come to a stop or not
Expand Down
92 changes: 40 additions & 52 deletions scripts/Reco/draw_spill.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,10 +29,15 @@
pad3.SetBottomMargin(0.15)
pad3.SetLeftMargin(0)
pad3.SetTopMargin(0.05)
pad4 = ROOT.TPad("p4", "p4", 0.0,0.,1,0.21)
pad4 = ROOT.TPad("p4", "p4", 0.0,0.,0.7,0.21)
pad4.SetBottomMargin(0.35)
pad4.SetRightMargin(0.025)
pad4.SetRightMargin(0.05)
pad4.SetTopMargin(0.05)
pad5 = ROOT.TPad("p5", "p5", 0.7,0.,1,0.21)
pad5.SetBottomMargin(0.35)
pad5.SetRightMargin(0.05)
pad5.SetLeftMargin(0.15)
pad5.SetTopMargin(0.05)

def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_filename, only_true_tms_muons = False):
if not os.path.exists(input_filename): raise ValueError(f"Cannot find input_filename {input_filename}")
Expand All @@ -53,6 +58,11 @@ def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_
bound_x_min = -4.0
bound_x_max = 4.0

slice_colors = [ROOT.kRed, ROOT.kBlue, ROOT.kGreen, ROOT.kOrange, ROOT.kPink, ROOT.kCyan, ROOT.kMagenta, ROOT.kAzure]

# TODO save this info in the file
spill_time = 1.2e9

font_size = 0.15

markers = []
Expand Down Expand Up @@ -88,17 +98,28 @@ def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_
# First loop through all the slices and draw one overall spill
for current_spill_number in range(max_n_spills):
htime = ROOT.TH1D(f"htime_{current_spill_number}", ";Time (ns);Total Hit E (MeV)", 500, -200, 12000)
htime.GetYaxis().SetTitleOffset(0.2)
htime.GetYaxis().SetTitleOffset(0.35)
htime.GetXaxis().SetTitleSize(font_size)
htime.GetXaxis().SetLabelSize(font_size)
htime.GetYaxis().SetTitleSize(font_size*0.8)
htime.GetYaxis().SetLabelSize(font_size*0.8)
htime.SetFillColor(ROOT.kBlack)
htime.SetLineColor(ROOT.kBlack)

max_hit_energy = 50
henergy = ROOT.TH1D(f"henergy_{current_spill_number}", ";Hit E (MeV);N Hits", 50, 0, max_hit_energy)
henergy.GetYaxis().SetTitleOffset(0.5)
henergy.GetXaxis().SetTitleSize(font_size)
henergy.GetXaxis().SetLabelSize(font_size)
henergy.GetYaxis().SetTitleSize(font_size*0.8)
henergy.GetYaxis().SetLabelSize(font_size*0.8)
henergy.SetFillColor(ROOT.kBlack)
henergy.SetLineColor(ROOT.kBlack)

total_energy = 0
time_high = float("-inf")
time_low = float("inf")
n_hits = 0

markers = []
text_to_draw = []
Expand Down Expand Up @@ -142,7 +163,7 @@ def inside_tms(x, y, z):
if not end_inside_tms: continue
print(f"Muon Start XYZ: ({mx:0.2f}, {my:0.2f}, {mz:0.2f})\tMuon End XYZ: ({mdx:0.2f}, {mdy:0.2f}, {mdz:0.2f})\tstart_inside_tms: {start_inside_tms}\tend_inside_tms: {end_inside_tms},\tPDG: {truth.LeptonPDG}")

print(f"Event {i} has {event.nHits} hits, {event.nClusters} clusters, and {event.nLines} lines.")
print(f"Event {i} has {event.nHits} hits, and {event.nLines3D} lines.")


for hit in range(event.nHits):
Expand All @@ -153,6 +174,9 @@ def inside_tms(x, y, z):
reco_hit_z = event.RecoHitPos[hit*4 + 2] / 1000.0
reco_hit_time = event.RecoHitPos[hit*4 + 3]
reco_hit_energy = event.RecoHitEnergy[hit]
reco_hit_slice = event.RecoHitSlice[hit]

n_hits += 1

time_high = max(reco_hit_time, time_high)
time_low = min(reco_hit_time, time_low)
Expand All @@ -161,64 +185,24 @@ def inside_tms(x, y, z):
e = int(min(255, 255 * reco_hit_energy / 10.0))
t = int(min(255, 255 * reco_hit_time / 10000.0))

htime.Fill(reco_hit_time, reco_hit_energy)
htime.Fill(reco_hit_time % spill_time, reco_hit_energy)
henergy.Fill(reco_hit_energy if reco_hit_energy < max_hit_energy else 0.95 * max_hit_energy)
x = reco_hit_z
y = reco_hit_x

marker = ROOT.TMarker(x, y, 21)
color = ROOT.TColor.GetColor(e, 32, 32)
#color = ROOT.TColor.GetColor(e, 32, 32)
color = ROOT.kGray if reco_hit_slice == 0 else slice_colors[reco_hit_slice % len(slice_colors)] + reco_hit_slice // len(slice_colors) - 3
marker.SetMarkerColor(color)
markers.append(marker)

if reco_hit_energy > 0:
total_energy += reco_hit_energy
elif reco_hit_energy < 0:
print(f"Found unexpected reco_hit_energy < 0: {reco_hit_energy}")

# Extract all hits in clusters and tracks that are not -999 for easier access later on
TotalHitsInClusters = sum(numpy.array([event.nHitsInCluster[i] for i in range(min(25, event.nClusters))]))
FilteredClusterHitPos = numpy.empty(TotalHitsInClusters*2)
j = 0
for i in range(len(event.ClusterHitPos)):
if event.ClusterHitPos[i] != -999:
FilteredClusterHitPos[j] = event.ClusterHitPos[i]
j += 1

usedClusters = 0
for cluster in range(min(25, event.nClusters)):
cluster_energy = event.ClusterEnergy[cluster]
cluster_time = event.ClusterTime[cluster]
cluster_pos_z = event.ClusterPosMean[cluster*2 + 0] / 1000.0
cluster_pos_x = event.ClusterPosMean[cluster*2 + 1] / 1000.0
cluster_pos_z_std_dev = event.ClusterPosStdDev[cluster*2 + 0] / 1000.0
cluster_pos_x_std_dev = event.ClusterPosStdDev[cluster*2 + 1] / 1000.0
#cluster_total_std_dev = math.sqrt(cluster_pos_z_std_dev**2 + cluster_pos_x_std_dev**2)
cluster_total_std_dev = max(cluster_pos_z_std_dev, cluster_pos_x_std_dev)
e = int(min(255, 255 * cluster_energy / 10.0))
t = int(min(255, 255 * cluster_time / 10000.0))
color = ROOT.kAzure #kBlack # ROOT.TColor.GetColor(e, t, 128)
x = cluster_pos_z
y = cluster_pos_x
#print("cluster:", cluster, e, t, 255 * cluster_energy / 10.0, 255 * cluster_time / 10000.0, color, x, y, cluster_total_std_dev)

marker = ROOT.TMarker(x, y, 53)
marker.SetMarkerColor(color)
marker.SetMarkerSize(20 * cluster_total_std_dev)
markers.append(marker)
marker = 0

for ClusterHit in range(event.nHitsInCluster[cluster]):
hit_z = FilteredClusterHitPos[usedClusters*2 + ClusterHit*2 + 0] / 1000.0
hit_x = FilteredClusterHitPos[usedClusters*2 + ClusterHit*2 + 1] / 1000.0

marker = ROOT.TMarker(hit_z, hit_x, 21)
marker.SetMarkerColor(ROOT.kAzure-8)
marker.SetMarkerSize(0.5)
markers.append(marker)
marker = 0
usedClusters += event.nHitsInCluster[cluster]

TotalHitsInTracks = sum(numpy.array([event.nHitsInTrack[i] for i in range(event.nLines)]))
'''TotalHitsInTracks = sum(numpy.array([event.nHitsInTrack[i] for i in range(event.nLines3D)]))
FilteredTrackHitPos = numpy.empty(TotalHitsInTracks*2)
j = 0
for i in range(len(event.TrackHitPos)):
Expand All @@ -227,7 +211,7 @@ def inside_tms(x, y, z):
j += 1
usedTracks = 0
for line in range(event.nLines):
for line in range(event.nLines3D):
#track_z = event.TrackHitPos[line*2 + 0] / 1000.0
#track_x = event.TrackHitPos[line*2 + 1] / 1000.0
track_z = event.FirstHoughHit[line*2 + 0] / 1000.0
Expand Down Expand Up @@ -277,12 +261,12 @@ def inside_tms(x, y, z):
marker.SetMarkerSize(0.5)
markers.append(marker)
marker = 0
usedTracks += event.nHitsInTrack[line]
usedTracks += event.nHitsInTrack[line]'''


minimum_energy_to_print = 0
if total_energy > minimum_energy_to_print:
print(f"Time range: {time_low:0.0f}-{time_high:0.0f} ns,\tEnergy: {total_energy:0.2f} MeV")
print(f"Spill {current_spill_number},\tTime range: {time_low:0.0f}-{time_high:0.0f} ns,\tEnergy: {total_energy:0.2f} MeV,\tN hits: {n_hits}")
pad1.cd()
haxis.Draw("colz axis")
for marker in markers:
Expand All @@ -292,10 +276,14 @@ def inside_tms(x, y, z):
# Draw the time slice
pad4.cd()
htime.Draw("hist")
# Draw the time slice
pad5.cd()
henergy.Draw("hist")

canvas.cd()
pad1.Draw()
pad4.Draw()
pad5.Draw()
output_filename = os.path.join(out_dir, f"{name}_{current_spill_number:03d}")
canvas.Print(output_filename + ".png")

Expand Down
Loading

0 comments on commit 14f8a6b

Please sign in to comment.