Skip to content

Commit

Permalink
Merge pull request #116 from DUNE/kleykamp_additional_truth_changes
Browse files Browse the repository at this point in the history
Kleykamp additional truth changes
  • Loading branch information
LiamOS authored Jul 11, 2024
2 parents cb3f4a7 + 1744bb6 commit 0f4bb23
Show file tree
Hide file tree
Showing 28 changed files with 3,120 additions and 299 deletions.
5 changes: 4 additions & 1 deletion app/CherryPickEvents.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -57,11 +57,14 @@ int main(int argc, char** argv) {
// Get the true neutrino vector from the gRooTracker object
int StdHepPdg[__EDEP_SIM_MAX_PART__];
double StdHepP4[__EDEP_SIM_MAX_PART__][4];
double StdHepX4[__EDEP_SIM_MAX_PART__][4];
gRoo->SetBranchStatus("*", false);
gRoo->SetBranchStatus("StdHepPdg", true);
gRoo->SetBranchStatus("StdHepP4", true);
gRoo->SetBranchStatus("StdHepX4", true);
gRoo->SetBranchAddress("StdHepPdg", StdHepPdg);
gRoo->SetBranchAddress("StdHepP4", StdHepP4);
gRoo->SetBranchAddress("StdHepX4", StdHepX4);

// Get the event
TG4Event *event = NULL;
Expand Down Expand Up @@ -183,7 +186,7 @@ int main(int argc, char** argv) {
}

TMS_Event tms_event = TMS_Event(*event, true);
tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4);
tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, StdHepX4);

int pdg = tms_event.GetNeutrinoPDG();
double enu = tms_event.GetNeutrinoP4().E();
Expand Down
26 changes: 20 additions & 6 deletions app/ConvertToTMSTree.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -43,17 +43,19 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) {
TG4Event *event = NULL;
events->SetBranchAddress("Event", &event);
// Get the true neutrino vector from the gRooTracker object
// Nice list of vars to avoid future confusion (internet search for genie groottracker)
// https://internal.dunescience.org/doxygen/read__t2k__rootracker_8C.html
int StdHepPdg[__EDEP_SIM_MAX_PART__];
double StdHepP4[__EDEP_SIM_MAX_PART__][4];
double StdHepX4[__EDEP_SIM_MAX_PART__][4];
double EvtVtx[__EDEP_SIM_MAX_PART__][4];
if (gRoo){
gRoo->SetBranchStatus("*", false);
gRoo->SetBranchStatus("StdHepPdg", true);
gRoo->SetBranchStatus("StdHepP4", true);
gRoo->SetBranchStatus("StdHepX4", true);
gRoo->SetBranchStatus("EvtVtx", true);
gRoo->SetBranchAddress("StdHepPdg", StdHepPdg);
gRoo->SetBranchAddress("StdHepP4", StdHepP4);
gRoo->SetBranchAddress("StdHepX4", StdHepX4);
gRoo->SetBranchAddress("EvtVtx", EvtVtx);
}
// The global manager
TMS_Manager::GetInstance().SetFileName(filename);
Expand All @@ -74,6 +76,7 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) {
Timer.Start();

int i = 0;
int truth_info_entry_number = 0;
//N_entries = 500;

// Do we overlay events
Expand All @@ -100,7 +103,7 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) {

// Fill up truth information from the GRooTracker object
if (gRoo){
tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4);
tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, EvtVtx);
}

// Keep filling up the vector and move on to the next event
Expand All @@ -127,6 +130,18 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) {

int nslices = TMS_TimeSlicer::GetSlicer().RunTimeSlicer(tms_event);

// Check if this is not pileup
if (event->Primaries.size() == 1 && tms_event.GetNVertices() == 1) {
// Fill the info of the one and only true vertex in the spill
auto primary_vertex = event->Primaries[0];
int interaction_number = primary_vertex.GetInteractionNumber();
gRoo->GetEntry(interaction_number);
tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, EvtVtx);
}

TMS_TreeWriter::GetWriter().FillSpill(tms_event, truth_info_entry_number, nslices);
truth_info_entry_number += nslices;

// Could save per spill info here

//std::cout<<"Ran time slicer"<<std::endl;
Expand Down Expand Up @@ -159,8 +174,7 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) {
auto primary_vertex = event->Primaries[primary_vertex_id];
int interaction_number = primary_vertex.GetInteractionNumber();
gRoo->GetEntry(interaction_number);
tms_event_slice.FillTruthFromGRooTracker(StdHepPdg, StdHepP4);
tms_event_slice.FillAdditionalTruthFromGRooTracker(StdHepX4);
tms_event_slice.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, EvtVtx);
// And the lepton info
int lepton_index = -1;
int current_index = 0;
Expand Down
5 changes: 4 additions & 1 deletion app/DrawEvents.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -56,12 +56,15 @@ int main(int argc, char** argv) {
// Get the true neutrino vector from the gRooTracker object
int StdHepPdg[__EDEP_SIM_MAX_PART__];
double StdHepP4[__EDEP_SIM_MAX_PART__][4];
double StdHepX4[__EDEP_SIM_MAX_PART__][4];
if (gRoo) {
gRoo->SetBranchStatus("*", false);
gRoo->SetBranchStatus("StdHepPdg", true);
gRoo->SetBranchStatus("StdHepP4", true);
gRoo->SetBranchStatus("StdHepX4", true);
gRoo->SetBranchAddress("StdHepPdg", StdHepPdg);
gRoo->SetBranchAddress("StdHepP4", StdHepP4);
gRoo->SetBranchStatus("StdHepX4", true);
}

// Get the event
Expand Down Expand Up @@ -187,7 +190,7 @@ int main(int argc, char** argv) {

TMS_Event tms_event = TMS_Event(*event, true);
if (gRoo){
tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4);
tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, StdHepX4);
}

int pdg = tms_event.GetNeutrinoPDG();
Expand Down
15 changes: 12 additions & 3 deletions app/Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,6 @@ LD_SHARED=g++
ROOT_LIBS = $(shell root-config --evelibs)
ROOT_INC = $(shell root-config --cflags)

GEANT4_LIBS = $(shell geant4-config --libs)
GEANT4_INC = $(shell geant4-config --prefix)/include/

ifndef CLHEP_INC # If not defined, normal build.
CLHEP_INC = $(shell clhep-config --prefix)/include
CLHEP_LIBS = $(shell clhep-config --libs)
Expand All @@ -18,8 +15,20 @@ else # github CI passes CLHEP_INC as an environment var, use this to set lib dir
endif
endif

# Check if geant4-config is available
GEANT4_CONFIG := $(shell command -v geant4-config 2>/dev/null)

# Set GEANT4 and EDEP variables based on the availability of geant4-config
ifneq ($(GEANT4_CONFIG),)
GEANT4_LIBS = $(shell geant4-config --libs)
GEANT4_INC = $(shell geant4-config --prefix)/include/
EDEP_LIBS = -L $(EDEPSIM_LIB) -ledepsim_io $(CLHEP_LIBS)
EDEP_INC = -I$(EDEPSIM_INC) -I$(GEANT4_INC) -I$(CLHEP_INC)
else
EDEP_LIBS = -L $(EDEPSIM_LIB) -ledepsim_io
# Need CLHEP for unit conversion, can probably just move this into a standalone header to remove dependency?
EDEP_INC = -I$(EDEPSIM_INC) -I$(CLHEP_INC)
endif

# The TMS includes
TMS_INC = -I../src
Expand Down
20 changes: 19 additions & 1 deletion config/TMS_Default_Config.toml
Original file line number Diff line number Diff line change
Expand Up @@ -102,12 +102,30 @@

[Recon.Calibration]
EnergyCalibration = 0.08856 # MeV / PE


[Fiducial]
[Fiducial.TMS.Start]
X = -3300.0 # mm
Y = -2850.0 # mm
Z = 11362.0 # mm
[Fiducial.TMS.End]
X = 3300.0 # mm
Y = 160.0 # mm
Z = 18314.0 # mm
[Fiducial.LAr.Start]
X = -3478.48 # mm
Y = -2166.71 # mm
Z = 4179.24 # mm
[Fiducial.LAr.End]
X = 3478.48 # mm
Y = 829.282 # mm
Z = 9135.88 # mm

# Options for processing and saving truth information
# LightWeight reduces run time signficantly because we only care about TMS objects (no LAr)
[Truth]
LightWeight = false
LArFiducialCut = 200.0 # mm

# Draw PDF of "event display". Slows down reco considerably
[Applications]
Expand Down
24 changes: 12 additions & 12 deletions scripts/Reco/make_hists.py
Original file line number Diff line number Diff line change
Expand Up @@ -574,18 +574,18 @@ def validate_then_run(args):
if nfiles == 0:
raise ValueError(f"Did not find any files in {indir}")
print(f"Found {nfiles} files in {indir}")
if inlist != "":
# In this case the user specified a text file with the full paths
with open(inlist) as f:
file_data = f.read()
files_to_use = file_data.splitlines()
nfiles = len(files_to_use)
if nfiles == 0:
raise ValueError(f"Did not find any files in {inlist}")
print(f"Found {nfiles} files in {inlist}")
if infile != "":
# In this case, the user specified exactly one file. Usually they'd hadd many files together.
files_to_use = [infile]
if inlist != "":
# In this case the user specified a text file with the full paths
with open(inlist) as f:
file_data = f.read()
files_to_use = file_data.splitlines()
nfiles = len(files_to_use)
if nfiles == 0:
raise ValueError(f"Did not find any files in {inlist}")
print(f"Found {nfiles} files in {inlist}")
if infile != "":
# In this case, the user specified exactly one file. Usually they'd hadd many files together.
files_to_use = [infile]

outdir = args.outdir
if outdir == "":
Expand Down
78 changes: 78 additions & 0 deletions scripts/TTreeSlicer/LArSlicer.C
Original file line number Diff line number Diff line change
@@ -0,0 +1,78 @@
#include <TFile.h>
#include <TTree.h>
#include <TBranch.h>
#include <TTreeReader.h>
#include <TTreeReaderValue.h>
#include <TTreeReaderArray.h>

void selectEvents(const char* inputFile, const char* outputFilename, int nmax) {
// Open the input file
TFile *file = TFile::Open(inputFile);
if (!file || file->IsZombie()) {
printf("Error: Unable to open input file %s\n", inputFile);
return;
}

// Get the input tree
TTree *inputTree = (TTree*)file->Get("Truth_Info");
if (!inputTree) {
printf("Error: Unable to retrieve input Truth_Info\n");
file->Close();
return;
}

// Create a new file
TFile *outputFile = new TFile(outputFilename, "RECREATE");
if (!outputFile || outputFile->IsZombie()) {
printf("Error: Unable to create output file %s\n", outputFilename);
file->Close();
return;
}

// Clone the input tree
TTree *outputTree = inputTree->CloneTree(0); // Create an empty tree with the same branches

// Define the cut condition
// Example: Select events where positionX, positionY, and positionZ are within a certain range
const double LAr_Start_Exact[] = {-3478.48, -2166.71, 4179.24};
const double LAr_End_Exact[] = {3478.48, 829.282, 9135.88};
Double_t minX = LAr_Start_Exact[0], maxX = LAr_End_Exact[0];
Double_t minY = LAr_Start_Exact[1], maxY = LAr_End_Exact[1];
Double_t minZ = LAr_Start_Exact[2], maxZ = LAr_End_Exact[2];

// Set up TTreeReader for input tree
//TTreeReader reader(inputTree);
//inputTree->SetBranchStatus("*", true);
//TTreeReaderArray<float> X4(reader, "LeptonX4");

float X4[4];
inputTree->SetBranchAddress("LeptonX4", &X4);

// Loop over events
long n = 0;
while (n < inputTree->GetEntries() && (nmax < 0 || n < nmax)) {
inputTree->GetEntry(n);
// Access the position components
Double_t posX = X4[0];
Double_t posY = X4[1];
Double_t posZ = X4[2];
// Apply the cut condition
if (posX >= minX && posX <= maxX &&
posY >= minY && posY <= maxY &&
posZ >= minZ && posZ <= maxZ) {
outputTree->Fill(); // Fill the output tree with selected events
}
n += 1;
}

// Write the output tree to the output file
outputTree->Write();

// Clean up
outputFile->Close();
file->Close();
}

void LArSlicer(const char* inputFile, const char* outputFile, int nmax = -1) {
selectEvents(inputFile, outputFile, nmax);
}
Loading

0 comments on commit 0f4bb23

Please sign in to comment.