Kleykamp additional truth changes

app/CherryPickEvents.cpp

@@ -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;
@@ -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();

app/ConvertToTMSTree.cpp

@@ -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);
@@ -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
@@ -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
@@ -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;
@@ -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;

app/DrawEvents.cpp

@@ -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
@@ -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();

app/Makefile

@@ -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)
@@ -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

config/TMS_Default_Config.toml

@@ -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]

scripts/Reco/make_hists.py

@@ -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 == "":

scripts/TTreeSlicer/LArSlicer.C

@@ -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);
+}