diff --git a/.gitignore b/.gitignore
index c4d3cce..08756f9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,9 +1,12 @@
-*.root
+.vscode
*.o
*.swp
-*.root
*.log
.depend
+example/analysis
+example/analysis_fullrec
+example/data
+example/*.root
leaf/DataModelRootDict*
leaf/HKAstroAnalysis.*
leaf/DataModel
diff --git a/README.md b/README.md
index 0355e3a..b7afc7f 100644
--- a/README.md
+++ b/README.md
@@ -1,66 +1,81 @@
-# LEAF
-Low Energy Algoritm Framework
+# :seedling: LEAF :seedling:
+
+**L**ow **E**nergy **A**lgorithm **F**ramework
+
+
+# DESCRIPTION
This algorithm is an alternative and simple LE fitter than can be used for HyperK and SuperK.
~~~~~~~~~~~~~~~~~~~~~~~~~
2020/02/09: LEAF was convert as a C++ class and can be included in your code.
+~~~~~~~~~~~~~~~~~~~~~~~~~
-New compilation method:
- source RunAtStart.sh
- ./SetupDataModel.sh
- cd leaf/
- make clean; make
-
-In order to use the class in your code look at example/
-~~~~~~~~~~~~~~~~~~~~~~~~~
+# DOWNLOAD AND PRE-REQUISITES
+
+To download this repository use :
+
+```
+git clone --single-branch --branch develop_lite https://github.com/hyperk/LEAF.git
+```
+
+Then make sure the following pre-requisites are installed.
+1. ROOT v5r34 or superior (not tested for older versions, but might work).
+2. WCSim version compatible with your ROOT version.
+3. BONSAI installation (although LEAF can work without it).
+
+# COMPATIBILITY
-# Pre-requisite to use the code:
-1. BONSAI installation.
-2. ROOT v5r34 or superior (not tested for older versions, but might work).
-3. HKAstroAnalysis class is private and can be downloaded by SK collaborators on sukap cluster.
+Versions listed below have been tested so far.
-# Compatibility tested so far:
1. WCSim-hybrid version: for geometries "HyperK", "HyperK_mPMT", "HyperK_HybridmPMT", "HyperK_HybridmPMT10PC"
2. In general, with all WCSim-hybrid geometries using whether BoxandLine20inchHQE or PMT3inchR14374 PMTs.
3. With official HK WCSim: Ask G. Pronost.
-# How to:
-1. Source RunAtStart.sh after you updated your ROOT directory.
-2. Use the script ./SetupDataModel.sh to define the DataModel (if you have hk-AstroAnalysis, you should setup the global variable)
-2. Enter the leaf/ repository and make clean;make
-3. Enter the example repository and make clean;make
-4. One example of how to run the code is set in example: test_example.sh
-5. inputs PDF, input from WCSim can be downloaded on sukap cluster. Please untar them in the LEAF repository.
-6. You can use shell scripts in shell/ in order to run the fitter or launch on batch.
-
-# Useful scripts in ./macros and ./shell
-You can compile with GNUMake like following in ./macros:
+
+# INSTALLATION
+
+1. Go to libWCSIM/ repository and make softlinks to your WCSim installation.
+```
+$ cd libWCSIM/
+$ ln -s /path/to/your/WCSim/include/ include
+$ ln -s /path/to/your/WCSim/src/ src
+$ ln -s /path/to/your/WCSim/lib/libWCSimRoot.so libWCSimRoot.so
+$ ln -s /path/to/your/WCSim/lib/libWCSimRoot.so.1.12.xx libWCSimRoot.so.1.12.xx
```
+2. Go to the cloned LEAF repository and source the script RunAtStart.sh after making sure you are sourcing the proper ROOT directory.
+```
+$ cd /path/to/the/cloned/repository/LEAF
+$ source ./RunAtStart.sh
+```
+3. Enter the leaf/ repository and `make clean; make`
+4. Enter the example/ repository and `make clean; make`
+5. One example of how to run the code is set in example: test_example.sh
+
+# TUNING FILES
+
+Tuning files are provided in the inputs/ repository (e.g. ./inputs/timePDF_Directionality_DRnew.root) and are made from 10 MeV electrons generated isotropically and uniformly in the tank.
+Alternatively, you can generate your own PDFs by following the commands below.
+
+1. Go to macros/ repository and compile AnalyzeWSHierarchy and ProducePDF with GNUMake.
+```
+$ cd ./macros
+$ make AnalyzeWSHierarchy
$ make ProducePDF
```
-## Making tuning file (e.g. ./inputs/timePDF_Directionality_DRnew.root)
-1. Produce plots by AnalyzeWSHierarchy: reads out WCSim output and makes plots.
-2. Produce time PDF (and angular PDF) by ProducePDF: uses plots made by AnalyzeWSHierarchy and generate PDFs for LEAF.
+2. Produce plots using AnalyzeWSHierarchy which reads out WCSim output. Optionally, you can index a range of events to read with -s (start index) and -e (end index). Use option -h if when using a hybrid HK geometry presenting both PMTs and mPMTs.
```
-$ AnalyzeWSHierarchy -f wcrim_hybrid.root -o plots.root
-$ ProducePDF -f plots.root -o PDF.root
+$ ./AnalyzeWSHierarchy -i wcrim.root -o plots.root -s 0 -e 1000 -h
```
-## Making generic plots
-- LEAFOutputAnalysisHybrid_leafclass: read LEAF output to produce generic plots. If one uses the master branch for LEAF, please use LEAFOutputAnalysisHybrid_master
+3. Produce time PDF and angular PDF using ProducePDF: uses plots made using AnalyzeWSHierarchy and generate PDFs for LEAF.
+```
+$ ./ProducePDF -i wcsim.root -o PDF.root -s 0 -e 1000 -v
+```
-## Shell scripts for analysis of large files
-You can refere shell scripts in ./shell in order to analyze many files.
-They are not working with latest LEAF class and its examples. They are just example how to analyze.
+Pathes to PDF files required as inputs for LEAF are hard-coded inside leaf/LeafSplines.cc, inside the LoadSplines() methods.
+Make sure these pathes are consistant with the files you generate on your own.
-- generateShellXX.c: this is a root macro which generates shell scripts to be submitted to sukap by LauncherXX.sh
- - generateShell.c
- - generateShell_analyzeWCSim.c
-- LauncherXX.sh: this submits jobs to sukap
- - Launcher.sh
- - Launcher_analyzeWCSim.sh
-- Merger_analyzeWCSim.sh: merge generated output by Launcher_analyzeWCSim.sh
diff --git a/SetupDataModel.sh b/SetupDataModel.sh
deleted file mode 100755
index 55ccbe6..0000000
--- a/SetupDataModel.sh
+++ /dev/null
@@ -1,16 +0,0 @@
-#!/bin/bash
-
-if [ ! -z ./leaf/DataModel ]; then
- echo "Remove previous link"
- rm ./leaf/DataModel
-fi
-
-if [ "$HK_ASTROANALYSIS_DIR" == "" ] || [ -z $HK_ASTROANALYSIS_DIR ]; then
- echo "HK_ASTROANALYSIS_DIR is not defined or doesn't exist ($HK_ASTROANALYSIS_DIR). Will use lite DataModel."
- ln -s ${PWD}/leaf/DataModel-lite ./leaf/DataModel
-else
- echo "HK_ASTROANALYSIS_DIR was found to be: $HK_ASTROANALYSIS_DIR its DataModel will be used"
- ln -s $HK_ASTROANALYSIS_DIR/DataModelRoot ./leaf/DataModel
-fi
-
-
diff --git a/example/GNUmakefile b/example/GNUmakefile
index ede518e..28a8b73 100644
--- a/example/GNUmakefile
+++ b/example/GNUmakefile
@@ -1,52 +1,62 @@
# Makefile by Guillaume Pronost for LEAF example @ 2020/02/09
OSNAME = $(shell uname -s)
+LEAFDIR = ..
include ../Makefile/Makefile.${OSNAME}
-
# set compiler options for ROOT
CXXFLAGS += $(shell root-config --cflags)
-CXXFLAGS += '-fPIC' -std=c++11 -Wall -Wpedantic -Wno-long-long
+CXXFLAGS += -fPIC -std=c++17 -Wall -Wpedantic -Wno-long-long
+ROOT_VERSION = $(shell root-config --version | cut -d'.' -f1,2)
+ifeq ($(shell echo "$(ROOT_VERSION) >= 6.0" | bc -l), 1)
+ HK_USE_ROOT7 = 1
+ CXXFLAGS += -DHK_USE_ROOT7
+endif
-INCFLAGS = -I. -I$(shell root-config --incdir)
-INCFLAGS += -I$(WCSIM_BUILD_DIR)/include/WCSim
-# INCFLAGS += -I$(WCSIMDIR)/include
-# INCFLAGS += -I$(BONSAIDIR)/bonsai
+INCFLAGS = -I.
+#INCFLAGS += -I$(shell root-config --incdir) #included in --cflags
+INCFLAGS += -I$(LEAFDIR)/libWCSIM/include
INCFLAGS += -I$(LEAFDIR)/leaf
-INCFLAGS += -I$(LEAFDIR)/leaf/DataModel
+INCFLAGS += -I$(LEAFDIR)/leaf/DataModel-lite
-# LIBS += -L${WCSIMDIR} -lWCSimRoot
-LIBS += -L${WCSIM_BUILD_DIR}/lib -lWCSimRoot
+LIBS += -L$(LEAFDIR)/libWCSIM -lWCSimRoot
LIBS += $(shell root-config --libs) -lMinuit
# LIBS_BONSAI += -L${BONSAIDIR} -lWCSimBonsai
-LIBS_LEAF += -L${LEAFDIR}/lib -lDataModelLite -lHKManager -lLEAF #-lHKAstroAnalysis
-
+LIBS_LEAF += -L${LEAFDIR}/lib -lLEAF -lHKManager -lDataModelLite
-OBJECT = analysis
+OBJECT = analysis analysis_fullrec
-CXXFLAGS += -std=c++17
+all: ROOT_CHECK $(OBJECT)
-all: $(OBJECT)
+ROOT_CHECK:
+ifeq ($(HK_USE_ROOT7),1)
+ @echo '<< compile with ROOT7 >>'
+else
+ @echo '<< compile with old ROOT >>'
+ @echo $(DATAMODEL_SOURCES)
+endif
analysis: analysis.o
- @echo '<< compiling bin analysis >>'
- @$(CXX) -g $(CXXFLAGS) -o $@ $^ $(LIBS) $(LIBS_BONSAI) $(LIBS_LEAF)
-
-
+ @echo '<< compiling bin $@ >>'
+ @$(CXX) -g $(CXXFLAGS) $^ $(LIBS) $(LIBS_LEAF) -o $@
+
+analysis_fullrec: analysis_fullrec.o
+ @echo '<< compiling bin $@ >>'
+ @$(CXX) -g $(CXXFLAGS) $^ $(LIBS) $(LIBS_LEAF) -o $@
+
# default rules
.cc.o:
@echo '<< compiling' $< '>>'
@$(CXX) $(CXXFLAGS) $(INCFLAGS) -c $<
-
-%.o: %.cc %.hh
+
+.cpp.o:
@echo '<< compiling' $< '>>'
- @$(CXX) $(CXXFLAGS) $(INCFLAGS) -c -o $@ $<
+ @$(CXX) $(CXXFLAGS) $(INCFLAGS) -c $<
-%.o: %.cpp
+%.o: %.cpp %.hpp
@echo '<< compiling' $< '>>'
@$(CXX) $(CXXFLAGS) $(INCFLAGS) -c -o $@ $<
-
%.o: %.C %.h
@echo '<< compiling' $< '>>'
diff --git a/example/GNUmakefile_nobonsai b/example/GNUmakefile_nobonsai
deleted file mode 100644
index 888ec1b..0000000
--- a/example/GNUmakefile_nobonsai
+++ /dev/null
@@ -1,56 +0,0 @@
-# Makefile by Guillaume Pronost for LEAF example @ 2020/02/09
-
-OSNAME = $(shell uname -s)
-
-include ../Makefile/Makefile.${OSNAME}
-
-
-# set compiler options for ROOT
-CXXFLAGS += $(shell root-config --cflags)
-CXXFLAGS += '-fPIC' -std=c++11 -Wall -Wpedantic -Wno-long-long
-
-INCFLAGS = -I. -I$(shell root-config --incdir)
-INCFLAGS += -I$(WCSIMDIR)/include
-INCFLAGS += -I$(LEAFDIR)/leaf
-INCFLAGS += -I$(LEAFDIR)/leaf/DataModel
-
-LIBS += -L${WCSIMDIR} -lWCSimRoot
-LIBS += $(shell root-config --libs) -lMinuit
-LIBS_LEAF += -L${LEAFDIR}/lib -lDataModelLite -lHKManager -lLEAF #-lHKAstroAnalysis
-
-
-OBJECT = analysis
-
-all: $(OBJECT)
-
-analysis: analysis.o
- @echo '<< compiling bin analysis >>'
- @$(CXX) -g $(CXXFLAGS) -o $@ $^ $(LIBS) $(LIBS_LEAF)
-
-
-# default rules
-.cc.o:
- @echo '<< compiling' $< '>>'
- @$(CXX) $(CXXFLAGS) $(INCFLAGS) -c $<
-
-%.o: %.cc %.hh
- @echo '<< compiling' $< '>>'
- @$(CXX) $(CXXFLAGS) $(INCFLAGS) -c -o $@ $<
-
-%.o: %.cpp
- @echo '<< compiling' $< '>>'
- @$(CXX) $(CXXFLAGS) $(INCFLAGS) -c -o $@ $<
-
-
-%.o: %.C %.h
- @echo '<< compiling' $< '>>'
- @$(CXX) $(CXXFLAGS) $(INCFLAGS) -c -o $@ $<
-
-# cleaner
-clean:
- @echo '<< cleaning >>'
- @rm -f ./*.o
- @rm -f ./lib*.so
- @rm -f ./bin/*
- @rm -f $(OBJECT)
-
diff --git a/example/analysis.cpp b/example/analysis.cpp
index b2db7ef..bd43f8d 100644
--- a/example/analysis.cpp
+++ b/example/analysis.cpp
@@ -4,341 +4,203 @@
/** Date: Febuary 10th 2020 **/
/** Desc: Example application code for Benjamin's Low-E Fitter for Hyper-K **/
/*********************************************************************************/
-#include "analysis.h"
+#include "analysis.hpp"
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///* Some Utility Functions
-
-double calculateToWall(double R, double h, const std::vector& position, const std::vector& direction)
-{
- double x = position[0], y = position[1], z = position[2];
- double dx = direction[0], dy = direction[1], dz = direction[2];
-
- double A = dx * dx + dy * dy;
- double B = 2 * (x * dx + y * dy);
- double C = x * x + y * y - R * R;
-
- double discriminant = B * B - 4 * A * C;
-
- double t_cylinder = -1;
- if (discriminant >= 0) {
- double t1 = (-B - sqrt(discriminant)) / (2 * A);
- double t2 = (-B + sqrt(discriminant)) / (2 * A);
-
- if (t1 > 0) t_cylinder = t1;
- if (t2 > 0 && (t2 < t1 || t_cylinder < 0)) t_cylinder = t2;
- }
-
- //Calculate the t values for intersection with top and bottom
- double t_top = (h / 2 - z) / dz;
- double t_bottom = (-h / 2 - z) / dz;
-
- //Check if the intersections are on the top or bottom
- double x_top = x + t_top * dx;
- double y_top = y + t_top * dy;
- double x_bottom = x + t_bottom * dx;
- double y_bottom = y + t_bottom * dy;
-
- if (x_top * x_top + y_top * y_top > R * R) t_top = -1; //Invalid if outside radius
- if (x_bottom * x_bottom + y_bottom * y_bottom > R * R) t_bottom = -1; //Invalid if outside radius
-
- //Find the smallest positive t that corresponds to a valid intersection (side or top/bottom)
- double t_wall = -1;
- if (t_cylinder > 0) t_wall = t_cylinder;
- if (t_top > 0 && (t_top < t_wall || t_wall < 0)) t_wall = t_top;
- if (t_bottom > 0 && (t_bottom < t_wall || t_wall < 0)) t_wall = t_bottom;
-
- return t_wall;
-}
-
-double calculateDWall(double R, double h, const std::vector& position)
-{
- double x = position[0], y = position[1], z = position[2];
- double distanceToAxis = sqrt(x * x + y * y);
- double dSide = fabs(distanceToAxis - R);
- double dTop = fabs(z - h / 2);
- double dBottom = fabs(z + h / 2);
- return std::min({dSide, dTop, dBottom});
-}
-
-double angleBetween3DVectors(const std::vector& a, const std::vector& b)
-{
- if(a.size() < 3 || b.size() < 3)
- {
- throw std::invalid_argument("Vectors must be of size at least 3.");
- }
- // double dotProduct = dot(a,b);
-
- // Magnitudes
- double magA = std::sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
- double magB = std::sqrt(b[0]*b[0] + b[1]*b[1] + b[2]*b[2]);
-
- if (magA == 0 || magB == 0) {
- throw std::invalid_argument("Zero-length vector.");
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+arguments FetchInput(int argc, char** argv) {
+ int c = -1;
+ arguments arglist;
+
+ //Input in c the argument (-f etc...) and in optarg the next argument
+ //When the above test becomes -1, it means it fails to find a new argument
+ std::cout << std::endl;
+ while( (c = getopt(argc, argv, "i:o:d:h:s:e:v")) != -1 )
+ {
+ switch(c)
+ {
+ //Input file name
+ case 'i':
+ arglist.inputFile = optarg;
+ std::cout << "Input WCSim file: " << arglist.inputFile << std::endl;
+ break;
+
+ //Output file name
+ case 'o':
+ arglist.outputFile = optarg;
+ if(arglist.outputFile == ""){ arglist.outputFile = "out.txt";}
+ std::cout << "Output root file: " << arglist.outputFile << std::endl;
+ break;
+
+ //Darknoise
+ case 'd':
+ arglist.darkNoise = atof(optarg);
+ std::cout << "Dark noise frequency: " << arglist.darkNoise << " kHz" << std::endl;
+ break;
+
+ //Darknoise hybrid
+ case 'h':
+ arglist.darkNoiseH = atof(optarg);
+ arglist.hybrid = true;
+ std::cout << "Dark noise frequency (hybrid geometry): " << arglist.darkNoiseH << " kHz" << std::endl;
+ break;
+
+ //Starting event
+ case 's':
+ arglist.startEvent = atoi(optarg);
+ if(arglist.startEvent < 0){arglist.startEvent = 0;}
+ std::cout << "Starting event #" << arglist.startEvent << std::endl;
+ break;
+
+ //Ending event
+ case 'e':
+ arglist.endEvent = atoi(optarg);
+ if(arglist.endEvent < 0){arglist.endEvent = 0;}
+ if(arglist.endEvent >= arglist.startEvent){std::cout << "Ending event #" << arglist.endEvent << std::endl;}
+ if(arglist.endEvent < arglist.startEvent){std::cout << "Ending event = last WCSim event" << std::endl;}
+ break;
+
+ //Warning
+ case 'v':
+ arglist.verbose = true;
+ std::cout << "VERBOSE option on" << std::endl;
+ break;
+ }
}
+ std::cout << std::endl;
- // Cosine of the angle
- double cosTheta = dot(a,b) / (magA * magB);
-
- // Clamp to [-1, 1] to avoid domain errors in acos
- cosTheta = std::max(-1.0, std::min(1.0, cosTheta));
-
- // Return angle in degrees
- return std::acos(cosTheta) * 180.0 / M_PI;
-}
-
-bool IsDarkRateHit(WCSimRootTrigger *fRootTrigger, WCSimRootCherenkovDigiHit *wcDigitHit)
-{
- bool isDR = false;
- std::vector rawhitphotonIDs = wcDigitHit->GetPhotonIds();
- if(rawhitphotonIDs.size() > 0)
- {
- double DR_Amount = 0;
- for(unsigned int i = 0; i < rawhitphotonIDs.size(); i++)
- {
- TObject *RawHitTimess;
- if (rawhitphotonIDs[i] >= 0 && rawhitphotonIDs[i] < (fRootTrigger->GetCherenkovHitTimes())->GetEntries())
- {
- RawHitTimess = (fRootTrigger->GetCherenkovHitTimes())->At(rawhitphotonIDs[i]);
- }
- else
- {
- std::cerr << "Error: Invalid photon ID index: " << rawhitphotonIDs[i] << std::endl;
- continue;
- }
- WCSimRootCherenkovHitTime *wcRawHitTimee = dynamic_cast(RawHitTimess);
- // std::cout << "Photon type : " << wcRawHitTimee->GetPhotonCreatorProcessName() << std::endl;
- if(wcRawHitTimee->GetPhotonCreatorProcessName() == "darkNoise") DR_Amount+=1.0;
- }
- DR_Amount = DR_Amount / rawhitphotonIDs.size();
- if(DR_Amount > 0.1) isDR = true; //? what threshold to consider a hit as darkRate ?
- }
- else
- {
- std::cout << "Error: No photon IDs found for this hit : " << wcDigitHit->GetTubeId() << std::endl;
- }
- return isDR;
-}
-
-double GetResidualTime(const std::vector& origin, double originTime, WCSimRootCherenkovDigiHit *wcDigitHit, TimeDelta fTimeCorrection, double fLfTriggerTime)
-{
- float fCVacuum = 3e8 * 1e2 / 1e9; // speed of light, in centimeter per ns.
- float fNIndex = 1.373; // 1.385;//1.373;//refraction index of water
- double fLightSpeed = fCVacuum / fNIndex;
- double HitT = wcDigitHit->GetT() + fLfTriggerTime;
- WCSimRootPMT pmt = fLeafGeometry->GetPMT(wcDigitHit->GetTubeId() - 1, false);
- double PMTpos[3];
- for (int j = 0; j < 3; j++) PMTpos[j] = pmt.GetPosition(j);
- double distance = Astro_GetDistance(PMTpos, origin);
- double tof = distance / fLightSpeed;
- double hitTime = (fTimeCorrection + HitT) / TimeDelta::ns;
- return hitTime - tof - originTime;
-}
-
-double GetDistanceToNeighbors(WCSimRootTrigger *fRootTrigger, WCSimRootCherenkovDigiHit *wcDigitHit, int N_Neighbors)
-{
- std::vector distances;
-
- WCSimRootPMT pmt;
- pmt = fLeafGeometry->GetPMT(wcDigitHit->GetTubeId() - 1, false);
-
- std::vector PMTpos(3);
- for (int j = 0; j < 3; j++) PMTpos[j] = pmt.GetPosition(j);
-
- //loop over each hit, check if it is not the same as the one we are looking at
- for (int i = 0; i < fRootTrigger->GetNcherenkovdigihits(); i++)
- {
- WCSimRootCherenkovDigiHit *neighborHit = dynamic_cast((fRootTrigger->GetCherenkovDigiHits())->At(i));
-
- if(wcDigitHit->GetTubeId() == neighborHit->GetTubeId())
- continue;
-
- WCSimRootPMT neighborPMT;
- neighborPMT = fLeafGeometry->GetPMT(neighborHit->GetTubeId() - 1, false);
-
- std::vector neighborPMTpos(3);
- for (int j = 0; j < 3; j++) neighborPMTpos[j] = neighborPMT.GetPosition(j);
-
- double distance = calculateDistance(PMTpos, neighborPMTpos);
- distances.push_back(distance);
- }
-
- //* Sort the distances, keep the N_Neighbors closest ones
- std::sort(distances.begin(), distances.end());
- // size = N_Neighbors;
- if (distances.size() < (long unsigned int)N_Neighbors )
- {
- std::cerr << "Not enough neighbors found!" << std::endl;
- return -1;
- }
- double sum = 0;
- for (int i = 0; i < N_Neighbors; i++) sum += distances[i];
- return sum / N_Neighbors;
+ return arglist;
}
-///* MAIN EXECUTION
-
-int main(int argc, char **argv)
-{
- std::string sInputFile = "";
- std::string sOutputFile = "";
- int iNeededArgc = 3;
- double dDarkNoise = 0.; // kHz
- double dDarkNoiseHybrid = 0.; // kHz
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+/* MAIN EXECUTION */
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- iNeededArgc += 1;
+int main(int argc, char **argv) {
+ //Get arguments
+ arguments arglist = FetchInput(argc, argv);
+ std::string sInputFile = arglist.inputFile;
+ std::string sOutputFile = arglist.outputFile;
+ double dDarkNoise = arglist.darkNoise;
+ double dDarkNoiseHybrid = 0.;
#ifdef mPMT
- iNeededArgc += 1;
+ dDarkNoiseHybrid = arglist.darkNoiseH;
#endif
-
- if (argc == iNeededArgc)
- {
- int iArg = 1;
- sInputFile = argv[iArg];
- iArg += 1;
- sOutputFile = argv[iArg];
- iArg += 1;
- dDarkNoise = atof(argv[iArg]);
- iArg += 1;
-#ifdef mPMT
- dDarkNoiseHybrid = atof(argv[iArg]);
- iArg += 1;
-#endif
- }
- else
- {
-
- std::cout << "Synthax: " << argv[0] << " input output";
- std::cout << " DN_in_kHz_B&L";
-
-#ifdef mPMT
- std::cout << " DN_in_kHz_mPMT";
-#endif
- std::cout << std::endl;
- return 0;
- }
+ int firstEvents = arglist.startEvent; //We start counting at 0
+ int lastEvents = arglist.endEvent;
+ bool verbose = arglist.verbose;
// Read WCSim output
TFile *fInputFile = new TFile(sInputFile.c_str(), "READ");
// Get TTrees
- TTree *fInputTree = (TTree *)fInputFile->Get("wcsimT");
+ TTree *fInputTree = (TTree *)fInputFile->Get("wcsimT");
TTree *fInputGeoTree = (TTree *)fInputFile->Get("wcsimGeoT");
fLeafGeometry = 0;
fInputGeoTree->SetBranchAddress("wcsimrootgeom", &fLeafGeometry);
+ // Set Branch for PMTs
WCSimRootEvent *fIDevent = new WCSimRootEvent();
- // Set Branch
fInputTree->SetBranchAddress("wcsimrootevent", &fIDevent);
- // Set autodelete to avoid memory leak
- fInputTree->GetBranch("wcsimrootevent")->SetAutoDelete(kTRUE);
+ fInputTree->GetBranch("wcsimrootevent")->SetAutoDelete(kTRUE); // Set autodelete to avoid memory leak
- // WCSimRootEvent *fHybridevent = new WCSimRootEvent();
+ // Set Branch for mPMTs
#ifdef mPMT
- // Set Branche
- // fInputTree->SetBranchAddress("wcsimrootevent2", &fHybridevent);
- // Set autodelete to avoid memory leak
- fInputTree->GetBranch("wcsimrootevent2")->SetAutoDelete(kTRUE);
+ WCSimRootEvent *fHybridevent = new WCSimRootEvent();
+ fInputTree->SetBranchAddress("wcsimrootevent2", &fHybridevent);
+ fInputTree->GetBranch("wcsimrootevent2")->SetAutoDelete(kTRUE); // Set autodelete to avoid memory leak
#endif
+ // Set Branch for OD PMTs
#ifdef OD_ON
WCSimRootEvent *fODevent = new WCSimRootEvent();
- // Set Branche
fInputTree->SetBranchAddress("wcsimrootevent_OD", &fODevent);
- // Set autodelete to avoid memory leak
- fInputTree->GetBranch("wcsimrootevent_OD")->SetAutoDelete(kTRUE);
+ fInputTree->GetBranch("wcsimrootevent_OD")->SetAutoDelete(kTRUE); // Set autodelete to avoid memory leak
#endif
- // Read Geo
- fInputGeoTree->GetEntry(0);
-
// Create Output TTree
TFile *fOutputFile = new TFile(sOutputFile.c_str(), "RECREATE");
fOutputFile->SetCompressionLevel(2);
- TTree *fGeoTree = new TTree("wcsimGeoT", "Geometry TTree");
- TTree *fPrimaryTree = new TTree("Reduced", "Reduced TTree");
-
- // Set Branches
- SetGeoBranch(fGeoTree);
-
- FitterOutput leaf_output;
- SetCustomBranch(fPrimaryTree, leaf_output);
+ // Read Geo and set branch
+ fInputGeoTree->GetEntry(0);
+ TTree *fGeoTree = new TTree("wcsimGeoT", "Geometry TTree");
+ TTree *fPrimaryTree = new TTree("Reduced", "Reduced TTree");
+ fGeoTree->Branch("wcsimrootgeom", fLeafGeometry);
+ SetCustomBranch(fPrimaryTree);
fGeoTree->Fill();
// Get PMT Number:
int nPMT_ID = fLeafGeometry->GetWCNumPMT();
+ std::cout << " ID " << nPMT_ID << std::endl;
#ifdef OD_ON
int nPMT_OD = fLeafGeometry->GetODWCNumPMT();
+ std::cout << " OD " << nPMT_OD << std::endl;
#endif
int nMultPMT = fLeafGeometry->GetWCNumPMT(true);
-
- std::cout << " ID " << nPMT_ID << std::endl;
std::cout << " mPMT " << nMultPMT << std::endl;
+ // Initialize HK Manager geometry
HKManager::GetME()->SetGeometry(fLeafGeometry, dDarkNoise * 1e3, dDarkNoiseHybrid * 1e3);
// Initialize LEAF
- LEAF::GetME()->Initialize(HKManager::GetME()->GetGeometry()); // Loads the geometry in leaf and loads the pdfs
+ LEAF::GetME()->Initialize( HKManager::GetME()->GetGeometry(),
+ HKManager::GetME()->GetDarkNoise(),
+ HKManager::GetME()->GetGeometryPMT_ID(),
+ HKManager::GetME()->GetGeometryPMT_mPMT()); // Loads the geometry in leaf and loads the pdfs
//* for each config variable, check if it has been set, take the default variable if not
- int maxEvents = -1; // in cm, the step size for coarse grid search
- const char* envValue = std::getenv("nbOfEvents");
- if (envValue != nullptr)
- {
- maxEvents = std::stoi(envValue);
- std::cout << "max events is set to " << maxEvents << std::endl;
- } else std::cout << "Environment variable max events is not set. Using all events" << std::endl;
-
- envValue = std::getenv("maxHitsAngle");
- if (envValue != nullptr)
- {
+ const char* envValue = std::getenv("maxHitsAngle");
+ if (envValue != nullptr) {
maxHitAngle = std::stof(envValue);
std::cout << "max hit angle : " << maxHitAngle << "°" << std::endl;
- } else std::cout << "Environment variable max hit angle is not set. Using default value : " << maxHitAngle << std::endl;
+ }
+ else std::cout << "Environment variable max hit angle is not set. Using default value : " << maxHitAngle << std::endl;
envValue = std::getenv("N_Neighbors");
- if (envValue != nullptr)
- {
+ if (envValue != nullptr) {
N_Neighbors = std::stof(envValue);
std::cout << "N_Neighbors : " << N_Neighbors << std::endl;
- } else std::cout << "Environment variable N_Neighbors is not set. Using default value : " << N_Neighbors << std::endl;
+ }
+ else std::cout << "Environment variable N_Neighbors is not set. Using default value : " << N_Neighbors << std::endl;
envValue = std::getenv("maxDistanceToNeighbors");
- if (envValue != nullptr)
- {
+ if (envValue != nullptr) {
maxDistanceToNeighbors = std::stof(envValue);
std::cout << "maxDistanceToNeighbors : " << maxDistanceToNeighbors << "cm" << std::endl;
- } else std::cout << "Environment variable maxDistanceToNeighbors is not set. Using default value : " << maxDistanceToNeighbors << std::endl;
+ }
+ else std::cout << "Environment variable maxDistanceToNeighbors is not set. Using default value : " << maxDistanceToNeighbors << std::endl;
// Read Input Tree
- int nPrimaryEvents = fInputTree->GetEntries();
-
- if(maxEvents > 0) nPrimaryEvents = std::min(nPrimaryEvents, maxEvents);
-
int iWrite = 0;
int failAmount = 0;
+ int nPrimaryEvents = fInputTree->GetEntries();
+ if(lastEvents > 0) nPrimaryEvents = std::min(nPrimaryEvents, lastEvents);
TStopwatch timer;
timer.Reset();
timer.Start();
// Loop on Primary events
- for (int i = 0; i < nPrimaryEvents; i++)
- {
- // Reset Hit vector
+ for (int i=firstEvents ; i 10 ) break;
+ // Reset Hit vector and clock
HKManager::GetME()->ResetHitInfo();
// HKManager::GetME()->ResetSecondaryHitInfo();
-
- // if ( i%1000==0 ) {
timer.Stop();
- std::cout << "Event # = " << i << " / " << nPrimaryEvents << " ( " << timer.RealTime() << " )\n";
timer.Reset();
- timer.Start();
- //}
+ if(verbose) {
+ std::cout << "\n===========================================================================================================================================================" << std::endl;
+ std::cout << "Event # = " << i+1 << " / " << nPrimaryEvents << std::endl;
+ }
+
+ timer.Start();
fInputTree->GetEntry(i);
// Initialize output variables
@@ -348,37 +210,25 @@ int main(int argc, char **argv)
rawhit_num_noDN = 0;
true_particleId.clear();
- true_energy.clear();
- true_origin_X.clear();
- true_origin_Y.clear();
- true_origin_Z.clear();
- true_origin_T.clear();
+ true_energies.clear();
+ true_vertexes.clear();
+
leaf_Vertex.clear();
leaf_Dir.clear();
leaf_MyDir.clear();
leaf_QuickDir.clear();
-
+
digithit_pmtId.clear();
digithit_T.clear();
correctedDigithit_T.clear();
digithit_Q.clear();
- digithit_Angle.clear();
- digithit_NormAngle.clear();
- digithit_NeighborsDist.clear();
- hit_is_DR.clear();
Charge_PMT.clear();
- relativeAngle.clear();
- lf_relativeAngle.clear();
- hit_residual.clear();
rawhit_num = 0;
digithit_num = 0;
fLFTime = -9999;
- hit_5_15ns = 0;
- hit_50ns = 0;
-
Hit_ID = 0;
Hit_ID_50 = 0;
Hit_ID_200 = 0;
@@ -397,14 +247,6 @@ int main(int argc, char **argv)
fLastRawHit = 0;
- trueVertex = std::vector(4);
- trueDir = std::vector(3);
-
- lf_spatial_res = 0.;
- lf_Dir_res = 0.;
- lf_time_res = 0.;
- lf_energy_res = 0.;
-
bestTrigger = 0;
fLfTriggerTime = 0.;
@@ -412,18 +254,17 @@ int main(int argc, char **argv)
/* ID events */
/****************************************************************************************/
- fHit = 0;
- fHit_20 = 0;
- fHit_50 = 0;
+ fHit = 0;
+ fHit_20 = 0;
+ fHit_50 = 0;
fHit_200 = 0;
fHit_400 = 0;
// this is where we send all the informations of the event to leaf
// we return false if we don't want to process this event based on defined filters
- bool shouldProcess = AnalyseEvent(fIDevent, ID_EVENT);
+ bool shouldProcess = AnalyseEvent(fIDevent, PMTType::kID);
- if(!shouldProcess)
- {
+ if(!shouldProcess) {
failAmount++;
std::cout << "Event # = " << i << " cannot be processed " << std::endl;
continue;
@@ -436,26 +277,27 @@ int main(int argc, char **argv)
Hit_ID_400 = fHit_400;
/****************************************************************************************/
- /* mPMT events */
+ /* mPMT events */
/****************************************************************************************/
+#ifdef mPMT
fHit = 0;
fHit_20 = 0;
fHit_50 = 0;
fHit_200 = 0;
fHit_400 = 0;
-#ifdef mPMT
- // /*bool bmPMT =*/AnalyseEvent(fHybridevent, mPMT_EVENT);
-#endif
+ // /*bool bmPMT =*/AnalyseEvent(fHybridevent, PMTType::kmPMT);
+
Hit_mPMT = fHit;
Hit_mPMT_20 = fHit_20;
Hit_mPMT_50 = fHit_50;
Hit_mPMT_200 = fHit_200;
Hit_mPMT_400 = fHit_400;
+#endif
/****************************************************************************************/
- /* OD events */
+ /* OD events */
/****************************************************************************************/
#ifdef OD_ON
@@ -467,53 +309,73 @@ int main(int argc, char **argv)
// There should be a dedicated OD analyser, as many thing should be different than for ID
// Doesn't exist yet
- /*bool bOD =*/AnalyseODEvent(fODevent, OD_EVENT);
+ /*bool bOD =*/ AnalyseODEvent(fODevent, OD_EVENT);
Hit_OD = fHit;
Hit_OD_50 = fHit_50;
Hit_OD_200 = fHit_200;
Hit_OD_400 = fHit_400;
#endif
+
/****************************************************************************************/
- /* Benjamin Fitter */
+ /* Benjamin Quilain Fitter */
/****************************************************************************************/
TStopwatch timerLF;
timerLF.Reset();
timerLF.Start();
- // To be replaced by meaningful TriggerTime
- TimeDelta fDummyTrigger(0.);
+ LEAF::GetME()->LoadHitsCollection(HKManager::GetME()->GetHitCollection_ID(),HKManager::GetME()->GetHitCollection_mPMT());
+ fLeafOutput = LEAF::GetME()->MakeSequentialFit();
+ //std::cout << "(X, Y, Z, T): (" << fLeafOutput.vtx.X() << ", " << fLeafOutput.vtx.Y() << ", " << fLeafOutput.vtx.Z() << ", " << fLeafOutput.vtx.T() << ") NLL: " << fLeafOutput.vtx_nll << std::endl;
+ //std::cout << "(dX, dY, dZ, T): (" << fLeafOutput.dir.X() << ", " << fLeafOutput.dir.Y() << ", " << fLeafOutput.dir.Z() << ") NLL: " << fLeafOutput.dir_nll << std::endl;
- leaf_output = LEAF::GetME()->MakeSequentialFit(HKManager::GetME()->GetHitCollection(), fDummyTrigger);
// the fit method can be replaced by a joint fit (doesn't work better for now)
- // leaf_output = LEAF::GetME()->MakeJointFit(HKManager::GetME()->GetHitCollection(), fDummyTrigger);
- fOutputProps = LEAF::fOutputProps;
-
+ // fLeafOutput = LEAF::GetME()->MakeJointFit(HKManager::GetME()->GetHitCollection(), fDummyTrigger);
+ fLeafPerformances = LEAF::GetME()->GetPerformances();
+
timerLF.Stop();
fLFTime = timerLF.RealTime();
- std::cout << " LEAF took: " << timerLF.RealTime() << " for " << HKManager::GetME()->GetHitCollection()->Size() << " Hits";
- std::cout << " (Vtx Search : " << fOutputProps.Vtx_Search_ComputeTime << " , Vtx Minimize : " << fOutputProps.Vtx_Minimize_ComputeTime << ", Dir Search : " << fOutputProps.Dir_Search_ComputeTime << " , Dir Minimize : " << fOutputProps.Dir_Minimize_ComputeTime << " , Energy Fit : " << fOutputProps.Energy_Fit_ComputeTime << ")" << std::endl;
- bool validEvent = PostLeafAnalysis(fIDevent,ID_EVENT, leaf_output);
+ if (verbose) {
+ std::cout << " LEAF took: " << timerLF.RealTime() << " sec for " << HKManager::GetME()->GetHitCollection_ID()->size() + HKManager::GetME()->GetHitCollection_mPMT()->size() << " Hits";
+ std::cout << " (Vtx Search : " << fLeafPerformances.vtx_search_ct
+ << " , Vtx Minimize : " << fLeafPerformances.vtx_minimize_ct
+ << ", Dir Search : " << fLeafPerformances.dir_search_ct
+ << " , Dir Minimize : " << fLeafPerformances.dir_minimize_ct
+ << " , Energy Fit : " << fLeafPerformances.energy_fit_ct << ")" << std::endl;
+ std::cout << " True vertex: (" << true_vertex.X()
+ << ", " << true_vertex.Y()
+ << ", " << true_vertex.Z()
+ << ", " << true_vertex.T()<< ") [cm/cm/cm/ns] "
+ << " // direction: (" << true_dir.X()
+ << ", " << true_dir.Y()
+ << ", " << true_dir.Z()
+ << " // energy: " << true_energy << std::endl;
+ std::cout << " LEAF vertex: (" << fLeafOutput.vtx.X()
+ << ", " << fLeafOutput.vtx.Y()
+ << ", " << fLeafOutput.vtx.Z()
+ << ", " << fLeafOutput.vtx.T()<< ") [cm/cm/cm/ns]"
+ << " // direction: (" << fLeafOutput.dir.X()
+ << ", " << fLeafOutput.dir.Y()
+ << ", " << fLeafOutput.dir.Z()
+ << " // energy: " << fLeafOutput.energy << std::endl;
+ }
/****************************************************************************************/
/* Fill output tree */
/****************************************************************************************/
- if(validEvent)
- {
- fPrimaryTree->Fill();
- iWrite += 1;
- }
+
+ fPrimaryTree->Fill();
+ iWrite += 1;
}
- std::cout << "Event # = " << nPrimaryEvents << " / " << nPrimaryEvents << std::endl;
-
- std::cout << ((nPrimaryEvents - failAmount) / nPrimaryEvents) * 100 << "% of the events were processed" << std::endl;
+ std::cout << "\n===========================================================================================================================================================";
+ std::cout << "\n===========================================================================================================================================================\n" << std::endl;
+ std::cout << ((nPrimaryEvents - failAmount) / nPrimaryEvents) * 100 << "% of the events have been processed" << std::endl;
fOutputFile->cd();
-
fOutputFile->Write("", TObject::kOverwrite);
delete fPrimaryTree;
@@ -521,8 +383,9 @@ int main(int argc, char **argv)
return 1;
}
+
//* All the variables that will be kept in the output Tree
-void SetCustomBranch(TTree *fPrimaryTree, FitterOutput leaf_output)
+void SetCustomBranch(TTree *fPrimaryTree)
{
fPrimaryTree->Branch("eventId", &eventId, "eventId/I");
fPrimaryTree->Branch("nTrigger", &nTrigger, "nTrigger/I");
@@ -531,14 +394,9 @@ void SetCustomBranch(TTree *fPrimaryTree, FitterOutput leaf_output)
fPrimaryTree->Branch("digithit_num", &digithit_num, "digithit_num/I");
fPrimaryTree->Branch("true_particleId", &true_particleId);
- fPrimaryTree->Branch("true_origin_X", &true_origin_X);
- fPrimaryTree->Branch("true_origin_Y", &true_origin_Y);
- fPrimaryTree->Branch("true_origin_Z", &true_origin_Z);
- fPrimaryTree->Branch("true_origin_T", &true_origin_T);
- fPrimaryTree->Branch("true_energy", &true_energy);
- fPrimaryTree->Branch("true_dir", &trueDir, "true_dir[3]/D");
- fPrimaryTree->Branch("DWall", &dWall, "DWall/D");
- fPrimaryTree->Branch("ToWall", &toWall, "lf_wall/D");
+ fPrimaryTree->Branch("true_vertexes", &true_vertexes);
+ fPrimaryTree->Branch("true_energies", &true_energies);
+ fPrimaryTree->Branch("true_dir", &true_dir);
fPrimaryTree->Branch("ID_hits", &Hit_ID, "Hit_ID/I");
fPrimaryTree->Branch("ID_hits_50", &Hit_ID_50, "Hit_ID_50/I");
@@ -551,60 +409,33 @@ void SetCustomBranch(TTree *fPrimaryTree, FitterOutput leaf_output)
fPrimaryTree->Branch("mPMT_hits_200", &Hit_mPMT_200, "Hit_mPMT_200/I");
fPrimaryTree->Branch("mPMT_hits_400", &Hit_mPMT_400, "Hit_mPMT_400/I");
- fPrimaryTree->Branch("hit_5_15ns", &hit_5_15ns, "hit_5_15ns/I");
- fPrimaryTree->Branch("hit_50ns", &hit_50ns, "hit_50ns/I");
-
- fPrimaryTree->Branch("lf_vertex", &leaf_Vertex);
- fPrimaryTree->Branch("lf_NLL", &leaf_output.NLL, "lf_NLL/D");
- fPrimaryTree->Branch("lf_good", &leaf_output.NLLR, "lf_good/D");
- fPrimaryTree->Branch("lf_ctime", &fLFTime, "lf_ctime/D"); // Computation time
- fPrimaryTree->Branch("lf_energy", &leaf_output.Energy, "lf_energy/D");
- fPrimaryTree->Branch("lf_Dir", &leaf_Dir);
- fPrimaryTree->Branch("lf_MyDir", &leaf_MyDir);
- fPrimaryTree->Branch("lf_Quick_Dir", &leaf_QuickDir);
- fPrimaryTree->Branch("lf_Dir_NLL", &leaf_output.DNLL, "lf_Dir_NLL/D");
- fPrimaryTree->Branch("lf_DWall", &lf_dWall, "DWall/D");
- fPrimaryTree->Branch("lf_ToWall", &lf_ToWall, "lf_ToWall/D");
-
- fPrimaryTree->Branch("lf_spatial_res", &lf_spatial_res, "lf_spatial_res/D");
- fPrimaryTree->Branch("lf_time_res", &lf_time_res, "lf_time_res/D");
- fPrimaryTree->Branch("lf_Dir_res", &lf_Dir_res,"lf_Dir_res/D");
- fPrimaryTree->Branch("lf_MyDir_res", &lf_MyDir_res,"lf_MyDir_res/D");
- fPrimaryTree->Branch("lf_Quick_Dir_res", &lf_Quick_Dir_res,"lf_Quick_Dir_res/D");
- fPrimaryTree->Branch("lf_energy_res", &lf_energy_res, "lf_energy_res/D");
-
- fPrimaryTree->Branch("Leaf_ComputeTime", &fOutputProps.Leaf_ComputeTime, "Leaf_ComputeTime/D");
- fPrimaryTree->Branch("Vtx_Search_ComputeTime", &fOutputProps.Vtx_Search_ComputeTime, "Vtx_Search_ComputeTime/D");
- fPrimaryTree->Branch("Vtx_Minimize_ComputeTime", &fOutputProps.Vtx_Minimize_ComputeTime, "Vtx_Minimize_ComputeTime/D");
- fPrimaryTree->Branch("Dir_Quick_Search_ComputeTime", &fOutputProps.Dir_Quick_Search_ComputeTime, "Dir_Quick_Search_ComputeTime/D");
- fPrimaryTree->Branch("Dir_Search_ComputeTime", &fOutputProps.Dir_Search_ComputeTime, "Dir_Search_ComputeTime/D");
- fPrimaryTree->Branch("Dir_Minimize_ComputeTime", &fOutputProps.Dir_Minimize_ComputeTime, "Dir_Minimize_ComputeTime/D");
- fPrimaryTree->Branch("Energy_Fit_ComputeTime", &fOutputProps.Energy_Fit_ComputeTime, "Energy_Fit_ComputeTime/D");
-
+ fPrimaryTree->Branch("lf_vtx", &fLeafOutput.vtx);
+ fPrimaryTree->Branch("lf_vtx_nll", &fLeafOutput.vtx_nll);
+ fPrimaryTree->Branch("lf_dir", &fLeafOutput.dir);
+ fPrimaryTree->Branch("lf_dir_nll", &fLeafOutput.dir_nll);
+ fPrimaryTree->Branch("lf_energy", &fLeafOutput.energy);
+ fPrimaryTree->Branch("lf_nll", &fLeafOutput.nll_r);
+ fPrimaryTree->Branch("lf_ct", &fLeafPerformances.leaf_ct);
+ fPrimaryTree->Branch("lf_vtx_search_ct", &fLeafPerformances.vtx_search_ct);
+ fPrimaryTree->Branch("lf_vtx_minization_ct", &fLeafPerformances.vtx_minimize_ct);
+ fPrimaryTree->Branch("lf_dir_search_ct", &fLeafPerformances.dir_search_ct);
+ fPrimaryTree->Branch("lf_dir_minization_ct", &fLeafPerformances.dir_minimize_ct);
+
fPrimaryTree->Branch("rawTriggerTime", &rawTriggerTime, "rawTriggerTime/D");
- fPrimaryTree->Branch("TotalCharge", &leaf_output.TotalCharge, "TotalCharge/D");
- fPrimaryTree->Branch("RelativeAngle", &relativeAngle);
- fPrimaryTree->Branch("lf_RelativeAngle", &lf_relativeAngle);
+ fPrimaryTree->Branch("TotalCharge", &fLeafOutput.total_charge, "TotalCharge/D");
fPrimaryTree->Branch("DigiHitT", &digithit_T);
fPrimaryTree->Branch("CorrectedDigiHitT", &correctedDigithit_T);
fPrimaryTree->Branch("HitIsDR", &hit_is_DR);
fPrimaryTree->Branch("DigiHitQ", &digithit_Q);
- fPrimaryTree->Branch("DigiHitAngle", &digithit_Angle);
- fPrimaryTree->Branch("DigiHitNormAngle", &digithit_NormAngle);
- fPrimaryTree->Branch("DigiHitNeighborsDist", &digithit_NeighborsDist);
fPrimaryTree->Branch("DigiHitResidual", &hit_residual);
fPrimaryTree->Branch("Charge_PMT", &Charge_PMT);
}
-void SetGeoBranch(TTree *fGeoTree)
-{
- fGeoTree->Branch("wcsimrootgeom", fLeafGeometry);
-}
-bool AnalyseEvent(WCSimRootEvent *tEvent, int iEventType)
-{
- int iHybrid = 0;
+bool AnalyseEvent(WCSimRootEvent *tEvent, PMTType pmtType) {
+ const int iPMTType = static_cast(pmtType);
+
int nVertex = 0; // Number of Vertex in event
int nTrack = 0; // Number of Track in event
int nRawCherenkovHits = 0; // Number of Raw Cherenkov hits
@@ -615,19 +446,14 @@ bool AnalyseEvent(WCSimRootEvent *tEvent, int iEventType)
rawTriggerTime = 0.;
WCSimRootTrigger *fRootTrigger;
- TimeDelta fDummyTrigger(0.);
- TimeDelta fTimeCorrection = HKManager::GetME()->GetHitCollection()->timestamp - fDummyTrigger; //*still have no idea what this is
-
//* find best trigger (best number of hits)
int bestNumbOfHits = 0;
bestTrigger = -1;
- for (int iTrig = 0; iTrig < nTrigger; iTrig++)
- {
+ for (int iTrig = 0; iTrig < nTrigger; iTrig++) {
fRootTrigger = tEvent->GetTrigger(iTrig);
nDigitizedCherenkovHits = fRootTrigger->GetNcherenkovdigihits();
- if(nDigitizedCherenkovHits > bestNumbOfHits && fRootTrigger->GetTriggerInfo().size() >= 3)
- {
+ if(nDigitizedCherenkovHits > bestNumbOfHits && fRootTrigger->GetTriggerInfo().size() >= 3) {
bestNumbOfHits = nDigitizedCherenkovHits;
bestTrigger = iTrig;
}
@@ -639,8 +465,7 @@ bool AnalyseEvent(WCSimRootEvent *tEvent, int iEventType)
//* Retreive true infos (onlyh set by wcsim in the first trigger)
- for (int iTrig = 0; iTrig < tEvent->GetNumberOfEvents(); iTrig++)
- {
+ for (int iTrig = 0; iTrig < tEvent->GetNumberOfEvents(); iTrig++) {
fRootTrigger = tEvent->GetTrigger(iTrig);
#ifdef OLD_WCSIM
@@ -655,26 +480,18 @@ bool AnalyseEvent(WCSimRootEvent *tEvent, int iEventType)
nDigitizedCherenkovHits = fRootTrigger->GetNcherenkovdigihits();
//* Fist trigger is the one used to compute true Vertex & Dir
- if (iTrig == 0)
- {
+ if (iTrig == 0) {
// Loop on vertex
- for (int iVertex = 0; iVertex < nVertex; iVertex++)
- {
-
+ for (int iVertex = 0; iVertex < nVertex; iVertex++) {
#ifdef OLD_WCSIM
- for(int j=0; j<3; j++) trueVertex[j]=fRootTrigger->GetVtx(j);
- true_origin_X.push_back(fRootTrigger->GetVtx(0));
- true_origin_Y.push_back(fRootTrigger->GetVtx(1));
- true_origin_Z.push_back(fRootTrigger->GetVtx(2));
+ true_vertex = ROOT::Math::XYZTVector(fRootTrigger->GetVtx(0),fRootTrigger->GetVtx(1),fRootTrigger->GetVtx(2),fRootTrigger->GetVtx(3));
+ true_vertexes.push_back(true_vertex);
#else
- for(int j=0; j<3; j++) trueVertex[j]=fRootTrigger->GetVtxs(iVertex,j);
- true_origin_X.push_back(fRootTrigger->GetVtxs(iVertex, 0));
- true_origin_Y.push_back(fRootTrigger->GetVtxs(iVertex, 1));
- true_origin_Z.push_back(fRootTrigger->GetVtxs(iVertex, 2));
+ true_vertex = ROOT::Math::XYZTVector(fRootTrigger->GetVtxs(iVertex,0),fRootTrigger->GetVtxs(iVertex,1),fRootTrigger->GetVtxs(iVertex,2),fRootTrigger->GetVtxs(iVertex,3));
+ true_vertexes.push_back(true_vertex);
#endif
- if (iVertex * 2 > nTrack)
- {
+ if (iVertex * 2 > nTrack) {
std::cout << "ERROR: Vertex and Track number incompatible (nVertex: " << nVertex << " ; nTrack " << nTrack << ") " << std::endl;
continue;
}
@@ -683,27 +500,19 @@ bool AnalyseEvent(WCSimRootEvent *tEvent, int iEventType)
TObject *Track = (fRootTrigger->GetTracks())->At(iVertex * 2);
WCSimRootTrack *wcTrack = dynamic_cast(Track);
- if (wcTrack->GetParentId() == 0 && wcTrack->GetIpnu()==11)
- {
- for(int j=0; j<3; j++) trueDir[j] = wcTrack->GetPdir(j);
- Normalize(trueDir);
+ if (wcTrack->GetParentId() == 0 && wcTrack->GetIpnu()==11) {
+ true_dir = ROOT::Math::XYZVector(wcTrack->GetPdir(0),wcTrack->GetPdir(1),wcTrack->GetPdir(2)).Unit();
}
-
- toWall = calculateToWall(R, h, trueVertex, trueDir);
- dWall = calculateDWall(R, h, trueVertex);
-
+
+ true_energy = wcTrack->GetE();
true_particleId.push_back(wcTrack->GetIpnu());
- true_energy.push_back(wcTrack->GetE());
- true_origin_T.push_back(fRootTrigger->GetVtx(3));
+ true_energies.push_back(wcTrack->GetE());
}
// Send True Position to LEAF to make checks of success in the steps of leaf (should be removed for true data)
- std::vector vVtxTrue(3, 0.);
- vVtxTrue[0] = true_origin_X[0];
- vVtxTrue[1] = true_origin_Y[0];
- vVtxTrue[2] = true_origin_Z[0];
- SetTrueVertexInfo(vVtxTrue, true_origin_T[0]);
- SetTrueDirInfo(trueDir);
+
+ LEAF::GetME()->SetTrueVertexInfo(true_vertexes[0]);
+ LEAF::GetME()->SetTrueDirInfo(true_dir);
rawhit_num = nRawCherenkovHits;
}
@@ -715,120 +524,66 @@ bool AnalyseEvent(WCSimRootEvent *tEvent, int iEventType)
rawTriggerTime = triggerInfo[2];
startingCherenkovHitID = digithit_pmtId.size();
- std::vector times;
-
- //* check if trigger is too soon or too late compared to the event (trigger happened because of dark rate)
- double hitTimesStack = 0;
- for (int iDigitHit = 0; iDigitHit < nDigitizedCherenkovHits; iDigitHit++)
- {
- TObject *Hit = (fRootTrigger->GetCherenkovDigiHits())->At(iDigitHit);
- WCSimRootCherenkovDigiHit *wcDigitHit = dynamic_cast(Hit);
- double HitT = wcDigitHit->GetT() + fLfTriggerTime;
- hitTimesStack+=HitT;
- }
- double hitTmean = hitTimesStack/nDigitizedCherenkovHits;
- if(hitTmean > 2000.0 || hitTmean < -2000.0) return false; //* set these values to define the filter (currently filtering nothing)
-
- double AllQ = 0;
// Loop on Digitized Hits
- for (int iDigitHit = 0; iDigitHit < nDigitizedCherenkovHits; iDigitHit++)
- {
- //* Get and Compute Hit Informations
-
+ double AllQ = 0;
+ for (int iDigitHit = 0; iDigitHit < nDigitizedCherenkovHits; iDigitHit++) {
+ //* Get Hit Informations
TObject *Hit = (fRootTrigger->GetCherenkovDigiHits())->At(iDigitHit);
WCSimRootCherenkovDigiHit *wcDigitHit = dynamic_cast(Hit);
int pmtId = wcDigitHit->GetTubeId();
WCSimRootPMT pmt = fLeafGeometry->GetPMT(pmtId - 1, false);
- double HitT = wcDigitHit->GetT() + fLfTriggerTime;
- double HitQ = wcDigitHit->GetQ();
+ double HitT = wcDigitHit->GetT() + fLfTriggerTime;
+ double HitQ = wcDigitHit->GetQ();
int peForTube = wcDigitHit->GetQ();
- std::vector lfHitDir = std::vector(3);
- std::vector PMTOrientation = std::vector(3);
- std::vector vDir = std::vector(3);
- bool isDR = IsDarkRateHit(fRootTrigger, wcDigitHit);
- double hitResidual = GetResidualTime(trueVertex, true_origin_T[0], wcDigitHit, fTimeCorrection, fLfTriggerTime);
- for (int j = 0; j < 3; j++)
- {
- PMTOrientation[j] = pmt.GetOrientation(j);
- vDir[j] = pmt.GetPosition(j) - trueVertex[j];
- lfHitDir[j] = pmt.GetPosition(j);
- }
-
- Normalize(vDir);
- Normalize(lfHitDir);
- Normalize(PMTOrientation);
-
- std::vector neg_vDir = { -vDir[0], -vDir[1], -vDir[2] };
- double NormAngle = dot(neg_vDir, PMTOrientation);
-
- double hitRelativeAngle = dot(vDir, trueDir);
-
- //* Use that inforation to filter some hits
- //* we aim at killing dark rate hits
- //* ideally without using any information form the true event
-
- //* Direction Filter
- double angle = angleBetween3DVectors(vDir, trueDir);
- // if(angle > maxHitAngle) continue; // works best at 90°
-
- //* Filter base on distance to hit neighbors
- double distanceToNeighbor = GetDistanceToNeighbors(fRootTrigger, wcDigitHit, N_Neighbors);
- // if(distanceToNeighbor > maxDistanceToNeighbors) continue;
-
- //*DR FILTER
- if(isDR) continue; // ideal filter to check reconstruction without dark rate
- ///* fill infos to output root file
+ //* Fill output root file
AllQ += HitQ;
digithit_pmtId.push_back(pmtId);
digithit_T.push_back(HitT);
correctedDigithit_T.push_back(HitT - fLfTriggerTime - triggerInfo[1]);
digithit_Q.push_back(HitQ);
- hit_is_DR.push_back(isDR);
- digithit_Angle.push_back(angle);
- digithit_NormAngle.push_back((TMath::ACos(NormAngle))*180./TMath::Pi());
- digithit_NeighborsDist.push_back(distanceToNeighbor);
Charge_PMT.push_back(peForTube);
- relativeAngle.push_back((TMath::ACos(hitRelativeAngle))*180./TMath::Pi());
- digithit_Type.push_back(iEventType);
- hit_residual.push_back(hitResidual);
+ digithit_Type.push_back(iPMTType);
}
- for(long unsigned int j=0;j -5) hit_5_15ns++;
- if(hit_residual[j] < 50 && hit_residual[j] > -50) hit_50ns++;
- }
int iIdx_BS = 0;
nDigitizedCherenkovHits = digithit_pmtId.size();
digithit_num = nDigitizedCherenkovHits - startingCherenkovHitID;
+ std::vector times;
// Feed fitter
- for (int iDigitHit = startingCherenkovHitID; iDigitHit < nDigitizedCherenkovHits; iDigitHit++)
- {
- times.push_back(digithit_T[iDigitHit]);
-
- if (digithit_pmtId[iDigitHit] <= 0) std::cout << " Weird PMT ID " << digithit_pmtId[iDigitHit] << std::endl;
- HKManager::GetME()->AddHit(digithit_T[iDigitHit], digithit_Q[iDigitHit], iHybrid, digithit_pmtId[iDigitHit]);
-
- // Bonsai (do not store mPMT hits)
- if (iIdx_BS < 2000 && digithit_T[iDigitHit] < 4000. && digithit_T[iDigitHit] > -4000. && iHybrid == 0)
- {
- // std::cout << pmtId << " " << HitT << std::endl;
- bsCAB[iIdx_BS] = digithit_pmtId[iDigitHit];
- bsT[iIdx_BS] = digithit_T[iDigitHit]; // shift BS time is needed if interaction time is < 0, this needs to be considered
- bsQ[iIdx_BS] = digithit_Q[iDigitHit];
-
- iIdx_BS += 1;
+ if (pmtType == PMTType::kID) {
+ for (int iDigitHit = startingCherenkovHitID; iDigitHit < nDigitizedCherenkovHits; iDigitHit++) {
+ if (digithit_pmtId[iDigitHit] <= 0) std::cout << " Weird PMT ID " << digithit_pmtId[iDigitHit] << std::endl;
+
+ times.push_back(digithit_T[iDigitHit]);
+ HKManager::GetME()->AddHit_ID(digithit_T[iDigitHit], digithit_Q[iDigitHit], digithit_pmtId[iDigitHit]);
+
+ // Bonsai (do not store mPMT hits)
+ if (iIdx_BS < 2000 && digithit_T[iDigitHit] < 4000. && digithit_T[iDigitHit] > -4000.) {
+ // std::cout << pmtId << " " << HitT << std::endl;
+ bsCAB[iIdx_BS] = digithit_pmtId[iDigitHit];
+ bsT[iIdx_BS] = digithit_T[iDigitHit]; // shift BS time is needed if interaction time is < 0, this needs to be considered
+ bsQ[iIdx_BS] = digithit_Q[iDigitHit];
+
+ iIdx_BS += 1;
+ bsnhit[0] = iIdx_BS;
+ }
}
}
+ else {
+ for (int iDigitHit = startingCherenkovHitID; iDigitHit < nDigitizedCherenkovHits; iDigitHit++) {
+ if (digithit_pmtId[iDigitHit] <= 0) std::cout << " Weird PMT ID " << digithit_pmtId[iDigitHit] << std::endl;
+ times.push_back(digithit_T[iDigitHit]);
+ HKManager::GetME()->AddHit_mPMT(digithit_T[iDigitHit], digithit_Q[iDigitHit], digithit_pmtId[iDigitHit]);
+ }
+
+ }
- // Bonsai
- bsnhit[0] = iIdx_BS;
fHit = digithit_pmtId.size();
// Compute hits
@@ -838,8 +593,7 @@ bool AnalyseEvent(WCSimRootEvent *tEvent, int iEventType)
std::vector Hit_time_400;
std::sort(times.begin(), times.end());
- for (int iDigitHit = 0; iDigitHit < nDigitizedCherenkovHits; iDigitHit++)
- {
+ for (int iDigitHit = 0; iDigitHit < nDigitizedCherenkovHits; iDigitHit++) {
double HitT = times[iDigitHit];
Hit_time_20.push_back(HitT);
@@ -853,7 +607,7 @@ bool AnalyseEvent(WCSimRootEvent *tEvent, int iEventType)
// Count hit in 50 ns window
while (HitT - Hit_time_50[0] > 50.) Hit_time_50.erase(Hit_time_50.begin());
- // Count hit in 200 ns windowAnalyseEvent
+ // Count hit in 200 ns window
while (HitT - Hit_time_200[0] > 200.) Hit_time_200.erase(Hit_time_200.begin());
// Count hit in 400 ns window
@@ -865,49 +619,6 @@ bool AnalyseEvent(WCSimRootEvent *tEvent, int iEventType)
if ((unsigned int)fHit_400 < Hit_time_400.size()) fHit_400 = Hit_time_400.size();
}
}
- return true;
-}
-bool PostLeafAnalysis(WCSimRootEvent * tEvent, int iEventType, FitterOutput leaf_output)
-{
- if (!(true_origin_X.size() > 0 && true_origin_Y.size() > 0 && true_origin_Z.size() > 0)) return false;
-
- for(int j=0;j<4;j++) leaf_Vertex.push_back(leaf_output.Vtx[j]);
- for(int j=0;j<3;j++) leaf_Dir.push_back(leaf_output.Dir[j]);
- for(int j=0;j<3;j++) leaf_MyDir.push_back(leaf_output.MyDir[j]);
- for(int j=0;j<3;j++) leaf_QuickDir.push_back(leaf_output.Quick_Dir[j]);
-
- //* RESOLUTIONS COMPUTATION
- std::vector true_vertex = {true_origin_X[0], true_origin_Y[0], true_origin_Z[0]};
- lf_dWall = calculateDWall(R, h, leaf_output.Vtx);
- lf_ToWall = calculateToWall(R, h,leaf_output.Vtx, leaf_output.Dir);
- lf_spatial_res = calculateDistance(true_vertex, leaf_output.Vtx);
- lf_time_res = abs(leaf_output.Vtx[3]-true_origin_T[0]);
-
- if(leaf_output.Dir.size() == 3) lf_Dir_res = TMath::ACos(dot(leaf_output.Dir, trueDir))*180./TMath::Pi();
- if(leaf_output.MyDir.size() == 3) lf_MyDir_res = TMath::ACos(dot(leaf_output.MyDir, trueDir))*180./TMath::Pi();
- if(leaf_output.Quick_Dir.size() == 3) lf_Quick_Dir_res = TMath::ACos(dot(leaf_output.Quick_Dir, trueDir))*180./TMath::Pi();
-
- lf_energy_res = abs(leaf_output.Energy - true_energy[0]);
-
- //* ANGLE ANALYSIS
- WCSimRootTrigger * fRootTrigger = tEvent->GetTrigger(bestTrigger);
-
- for(int iDigitHit = startingCherenkovHitID; (long unsigned int)iDigitHit < digithit_pmtId.size(); iDigitHit++)
- {
- TObject *Hit = (fRootTrigger->GetCherenkovDigiHits())->At(iDigitHit);
- WCSimRootCherenkovDigiHit *wcDigitHit = dynamic_cast(Hit);
- WCSimRootPMT pmt = fLeafGeometry->GetPMT(wcDigitHit->GetTubeId() - 1, false);
-
- std::vector lfHitDir = std::vector(3);
- for(int j=0;j<3;j++) lfHitDir[j] = pmt.GetPosition(j) - leaf_output.Vtx[j];
- Normalize(lfHitDir);
-
- std::vector trueHitDir = std::vector(3);
- for(int j=0;j<3;j++) trueHitDir[j] = pmt.GetPosition(j) - trueVertex[j];
- Normalize(trueHitDir);
- double hitlfRelativeAngle = dot(trueHitDir, leaf_output.Dir);
- lf_relativeAngle.push_back((TMath::ACos(hitlfRelativeAngle))*180./TMath::Pi());
- }
return true;
-}
\ No newline at end of file
+}
diff --git a/example/analysis.h b/example/analysis.hpp
similarity index 75%
rename from example/analysis.h
rename to example/analysis.hpp
index 6d22c78..250d7fb 100644
--- a/example/analysis.h
+++ b/example/analysis.hpp
@@ -24,17 +24,17 @@
#include "WCSimRootGeom.hh"
#include "WCSimEnumerations.hh"
-#include "LEAF.hh"
-#include "LeafDefinitions.hh"
-#include "HKManager.hh"
+#include "LEAF.hpp"
+#include "LEAFUtilities.hpp"
+#include "HKManager.hpp"
#define OLD_WCSIM // To be used if WCSim version is older than 1.8 (i.e. without multi vertex)
// #define mPMT // To be used if you are using mPMT
-#define ID_EVENT 1
-#define OD_EVENT 2
-// #define mPMT_EVENT 0
-#define UNDEFINED_EVENT 0
+#define ID_EVENT PMTType::kID
+#define OD_EVENT PMTType::kOD
+#define mPMT_EVENT PMTType::kmPMT
+#define UNDEFINED_EVENT PMTType::kUndefined
//-----------------------------------------------------------------------------------------//
@@ -48,14 +48,14 @@ int usedTriggerId; // trigger Id used for the fit
double fLfTriggerTime;
int bestTrigger;
double rawTriggerTime;
+
std::vector true_particleId; // True particle Id (PDGId)
-std::vector true_energy; // True energy (MeV)
-std::vector true_origin_X; // True origin vertex (cm)
-std::vector true_origin_Y; // True origin vertex (cm)
-std::vector true_origin_Z; // True origin vertex (cm)
-std::vector true_origin_T; // True origin vertex (ns)
-std::vector trueDir;
-std::vector trueVertex;
+std::vector true_energies; // True energy (MeV)
+std::vector true_vertexes; // True origin vertex (cm)
+ROOT::Math::XYZVector true_dir;
+ROOT::Math::XYZTVector true_vertex;
+double true_energy;
+
double dWall =0;
double toWall =0;
int digithit_num; // Number of digit hit
@@ -110,11 +110,11 @@ std::vector Charge_PMT;
std::vector hit_residual; // List of residual times
//* Computation Times
-FitterOutputProps fOutputProps;
+LEAFStructures::FitterPerformances fLeafPerformances;
+LEAFStructures::FitterOutput fLeafOutput;
double fLFTime;
//* Other hits infos
-
int Hit_ID;
int Hit_ID_20;
int Hit_ID_50;
@@ -153,10 +153,38 @@ float bsT[2000];
float bsQ[2000];
int bsnhit[1];
+
//-----------------------------------------------------------------------------------------//
-void SetCustomBranch(TTree *fPrimaryTree, FitterOutput leaf_output);
-void SetCustomBranchInput(TTree *fPrimaryTree);
-void SetGeoBranch(TTree *fGeoTree);
-bool AnalyseEvent(WCSimRootEvent *tEvent, int iEventType);
-bool PostLeafAnalysis(WCSimRootEvent * tEvent, int iEventType, FitterOutput leaf_output);
+struct arguments
+{
+ std::string inputFile = "";
+ std::string outputFile = "";
+ double darkNoise = 4.2; // Dark noise frequency in Hz
+ double darkNoiseH = 0.; // Dark noise frequency in Hz when using the hybrid geometry
+ double timeshift = 0.; // Shift hit time
+ int startEvent = 0; // First event to analyze
+ int endEvent = 0; // Last event to analyze
+ bool hybrid = false;
+ bool verbose = false;
+};
+
+struct FitterAnalysis
+{
+ double Wall;
+ double Good;
+ int n50[3];
+ double dir[3][3];
+ double dir_goodness[3];
+ double dirKS[3];
+};
+
+
+//-----------------------------------------------------------------------------------------//
+
+void SetCustomBranch(TTree *fPrimaryTree);
+bool AnalyseEvent(WCSimRootEvent *tEvent, PMTType pmtType);
+arguments FetchInput(int argc, char* argv[]);
+
+
+
diff --git a/example/analysis_fullrec.cpp b/example/analysis_fullrec.cpp
new file mode 100644
index 0000000..886a1f1
--- /dev/null
+++ b/example/analysis_fullrec.cpp
@@ -0,0 +1,1196 @@
+/*********************************************************************************/
+/** analysis_fullrec.cpp **/
+/** Author: Benjamin Quilain (benjamin.quilain@llr.in2p3.fr) **/
+/** Date: April 10th 2025 **/
+/** Desc: Perform the vertex, direction and energy reconstruction **/
+/*********************************************************************************/
+
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+
+#include "TFile.h"
+#include "TTree.h"
+#include "TBranch.h"
+#include "TObject.h"
+#include "TSystem.h"
+#include "TRandom3.h"
+#include "TStopwatch.h"
+
+#include "WCSimRootEvent.hh"
+#include "WCSimRootGeom.hh"
+#include "WCSimEnumerations.hh"
+
+#include "LEAF.hpp"
+#include "HKManager.hpp"
+//#define WITH_BONSAI
+#ifdef WITH_BONSAI
+#include "WCSimBonsai.hh"
+#endif
+
+//-----------------------------------------------------------------------------------------//
+
+ int eventId; // event Id
+ int triggerId; // trigger Id (with event Id)
+
+ double TotalQ;
+
+ ROOT::Math::XYZVector particleDir(0,0,0);
+ struct FitterAnalysis {
+ double Wall;
+ double Good;
+ int n50[3];
+ double dir[3][3];
+ double dir_goodness[3];
+
+ double dirKS[3];
+ };
+
+ void Normalize(double a[3]){
+ double l;
+ l=sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
+ for(int j=0; j<3; j++){
+ a[j]/=l;
+ }
+ }
+
+ //True information
+ std::vector true_particleId; // True particle Id (PDGId)
+ std::vector true_energies; // True energy (MeV)
+ std::vector true_origins; // True origin vertex (cm)
+ std::vector relativeAngle;
+ std::vector lf_relativeAngle;
+
+ ROOT::Math::XYZVector lf_Dir;
+ ROOT::Math::XYZVector lf_Dir_interp;
+ double lf_Dir_res;
+ WCSimRootGeom *geotree = 0;
+
+ //Raw hit
+ std::vector rawhit_pmtId; // List of pmtId
+ std::vector rawhit_Type; // List of pmt Type
+ std::vector rawhit_T; // List of hit times
+ std::vector rawhit_Q; // List of hit charge'
+ std::vector rawhit_dark; // List of DarkNoise flag
+ int rawhit_num; // Number of rawhit
+ int rawhit_num_noDN; // Number of rawhit without DarkNoise
+
+
+ std::vector rawhit_pmtId_tmp; // List of pmtId
+ std::vector rawhit_Type_tmp; // List of pmt Type
+ std::vector rawhit_T_tmp; // List of hit times
+
+ // In case of mPMT mixed with B&L we don't want to re apply the Digitization on the other type PMTs
+ int fLastRawHit;
+
+
+ //Digitized hit
+ std::vector digithit_pmtId; // List of pmtId
+ std::vector digithit_Type; // List of pmt Type
+ std::vector digithit_T; // List of hit times
+ std::vector digithit_Q; // List of hit charge'
+ std::vector Charge_PMT;
+ std::vector digithit_dark; // List of DarkNoise flag
+ int digithit_num; // Number of digit hit
+ int digithit_num_noDN; // Number of digit hit without DarkNoise
+
+ //Bonsai output
+ float bs_vertex[4]; // Bonsai reconstructed vertex (cm)
+ float bs_good[3]; // Bonsai goodness
+ float bs_energy;
+
+ double fBSTime;
+ double fLFTime;
+ double lf_spatial_res;
+ double lf_jitter;
+ LEAFStructures::FitterOutput leaf_output;
+ FitterAnalysis leaf_output_ana;
+ FitterAnalysis bs_output_ana;
+
+ int Hit_ID;
+ int Hit_ID_50;
+ int Hit_ID_200;
+ int Hit_ID_400;
+
+ int Hit_mPMT;
+ int Hit_mPMT_50;
+ int Hit_mPMT_200;
+ int Hit_mPMT_400;
+
+ int Hit_OD;
+ int Hit_OD_50;
+ int Hit_OD_200;
+ int Hit_OD_400;
+
+ int fHit = 0;
+ int fHit_50 = 0;
+ int fHit_200 = 0;
+ int fHit_400 = 0;
+ double dWall =0;
+ double lf_dWall =0;
+ double lf_ToWall =0;
+ WCSimRootGeom * fGeometry;
+
+ double R = 3242.96;
+ double h = 6701.41;
+
+ // Input variable for Bonsai
+ int bsCAB[2000];
+ float bsT[2000];
+ float bsQ[2000];
+ int bsnhit[1];
+
+//-----------------------------------------------------------------------------------------//
+
+
+void SetCustomBranch(TTree* fPrimaryTree);
+void SetCustomBranchInput(TTree* fPrimaryTree);
+void SetGeoBranch(TTree* fGeoTree);
+
+// Apply analysis on the event
+bool AnalyseEvent(WCSimRootEvent * tEvent, PMTType iEventType, int i);
+bool AnalyselfEvent(WCSimRootEvent * tEvent, PMTType iEventType);
+TH1D *hRelativeAngle = nullptr;
+TH1D *lf_hRelativeAngle = nullptr;
+TH1D *hRelativeAngleCos = nullptr;
+TH1D *lf_hRelativeAngleCos = nullptr;
+TH1D *ToWall_Charge = nullptr;
+TH2D *ToWall_Charge_2D = nullptr;
+TH1D* Sum_HitQ = nullptr;
+TH1D* Count_Hits = nullptr;
+
+TH2D* hToWall_RelAngle_Charge = new TH2D("ToWall_RelAngle_Charge", "ToWall vs Rel. Angle and Charge", 100, 0, 9000, 360, 0, 180);
+TProfile* pAngle_ToWall = new TProfile("pAngle_ToWall", "Angle vs ToWall", 360, 0, 180);
+
+TH2D* hdWall_TotalCharge = new TH2D("DWall_TotalCharge", "dWall vs Total Charge", 150, 0, 3500, 50, 0, 400);
+TH2D* hdWall_TotalCharge_corr = new TH2D("DWall_TotalChargeCorr", "dWall vs Corrected Total Charge", 150, 0, 3500, 50, 0, 400);
+TH1D* hdWall_TotalCharge_hist = new TH1D("DWall_TotalCharge_hist", "dWall vs Total Charge (Hist)", 150, 0, 3500);
+TH1D* Sum_HitQ_hdWall = new TH1D("Sum_HitQ", "Sum of HitQ per Wall bin", 150, 0, 3500);
+TH1D* Count_Hits_hdWall = new TH1D("Count_Hits", "Count of Hits per Wall bin", 150, 0, 3500);
+
+
+
+TH1D* hAngle_HitQ = new TH1D("Angle_HitQ", "Angle vs Hit Charge", 360, 0, 90);
+
+void initializeHistograms() {
+ hToWall_RelAngle_Charge->GetXaxis()->SetTitle("ToWall [cm]");
+ hToWall_RelAngle_Charge->GetYaxis()->SetTitle("Relative Angle [degrees]");
+ hToWall_RelAngle_Charge->GetZaxis()->SetTitle("Charge [p.e.]");
+
+ hdWall_TotalCharge->GetXaxis()->SetTitle("dWall [cm]");
+ hdWall_TotalCharge->GetYaxis()->SetTitle("Total Charge [p.e.]");
+
+ hdWall_TotalCharge_corr->GetXaxis()->SetTitle("dWall [cm]");
+ hdWall_TotalCharge_corr->GetYaxis()->SetTitle("Corrected Total Charge [p.e.]");
+
+ hdWall_TotalCharge_hist->GetXaxis()->SetTitle("dWall [cm]");
+ hdWall_TotalCharge_hist->GetYaxis()->SetTitle("Total Charge [p.e.]");
+
+ hAngle_HitQ->GetXaxis()->SetTitle("Angle [degrees]");
+ hAngle_HitQ->GetYaxis()->SetTitle("Hit Charge [p.e.]");
+}
+
+double calculateToWall(double R, double h, const ROOT::Math::XYZVector& position, const ROOT::Math::XYZVector& direction) {
+ double x = position.X(), y = position.Y(), z = position.Z();
+ double dx = direction.X(), dy = direction.Y(), dz = direction.Z();
+
+ double A = dx * dx + dy * dy;
+ double B = 2 * (x * dx + y * dy);
+ double C = x * x + y * y - R * R;
+
+ double discriminant = B * B - 4 * A * C;
+
+ double t_cylinder = -1;
+ if (discriminant >= 0) {
+ double t1 = (-B - sqrt(discriminant)) / (2 * A);
+ double t2 = (-B + sqrt(discriminant)) / (2 * A);
+
+ if (t1 > 0) t_cylinder = t1;
+ if (t2 > 0 && (t2 < t1 || t_cylinder < 0)) t_cylinder = t2;
+ }
+
+ //Calculate the t values for intersection with top and bottom
+ double t_top = (h / 2 - z) / dz;
+ double t_bottom = (-h / 2 - z) / dz;
+
+ //Check if the intersections are on the top or bottom
+ double x_top = x + t_top * dx;
+ double y_top = y + t_top * dy;
+ double x_bottom = x + t_bottom * dx;
+ double y_bottom = y + t_bottom * dy;
+
+ if (x_top * x_top + y_top * y_top > R * R) t_top = -1; //Invalid if outside radius
+ if (x_bottom * x_bottom + y_bottom * y_bottom > R * R) t_bottom = -1; //Invalid if outside radius
+
+ //Find the smallest positive t that corresponds to a valid intersection (side or top/bottom)
+ double t_wall = -1;
+ if (t_cylinder > 0) t_wall = t_cylinder;
+ if (t_top > 0 && (t_top < t_wall || t_wall < 0)) t_wall = t_top;
+ if (t_bottom > 0 && (t_bottom < t_wall || t_wall < 0)) t_wall = t_bottom;
+
+ return t_wall;
+}
+
+double calculateDWall(double R, double h, const ROOT::Math::XYZVector& position) {
+ double x = position.X(), y = position.Y(), z = position.Z();
+ double distanceToAxis = sqrt(x * x + y * y);
+ double dSide = fabs(distanceToAxis - R);
+ double dTop = fabs(z - h / 2);
+ double dBottom = fabs(z + h / 2);
+ return std::min({dSide, dTop, dBottom});
+}
+
+
+
+
+
+
+int main(int argc, char** argv){
+
+ initializeHistograms();
+
+
+
+
+ std::string sInputFile = "";
+ std::string sOutputFile = "";
+
+ int iNeededArgc = 5;
+ double dDarkNoise = 0.; //kHz
+ double dDarkNoiseHybrid = 0.; //kHz
+ iNeededArgc += 1;
+ int StartEvent;
+ int NumberEvents;
+#ifdef mPMT
+ iNeededArgc += 1;
+#endif
+
+ if ( argc == iNeededArgc ) {
+ int iArg = 1;
+ sInputFile = argv[iArg]; iArg += 1;
+ sOutputFile = argv[iArg]; iArg += 1;
+ dDarkNoise = atof(argv[iArg]); iArg += 1;
+#ifdef mPMT
+ dDarkNoiseHybrid = atof(argv[iArg]); iArg += 1;
+#endif
+ StartEvent = atoi(argv[iArg]); iArg += 1;
+ NumberEvents = atoi(argv[iArg]); iArg += 1;
+
+ }
+ else {
+
+ std::cout << "Synthax: " << argv[0] << " input output";
+ std::cout << " DN_in_kHz_B&L";
+
+#ifdef mPMT
+ std::cout << " DN_in_kHz_mPMT";
+#endif
+ std::cout << std::endl;
+ return 0;
+ }
+
+
+ // Read WCSim output
+ TFile* fInputFile = new TFile(sInputFile.c_str(),"READ");
+
+ // Get TTrees
+ TTree* fInputTree = (TTree*) fInputFile->Get("wcsimT");
+ TTree *fInputGeoTree = (TTree*) fInputFile->Get("wcsimGeoT");
+
+ fGeometry = 0;
+ fInputGeoTree->SetBranchAddress("wcsimrootgeom", &fGeometry);
+
+ WCSimRootEvent * fIDevent = new WCSimRootEvent();
+ // Set Branche
+ fInputTree->SetBranchAddress("wcsimrootevent" ,&fIDevent );
+ // Set autodelete to avoid memory leak
+ fInputTree->GetBranch("wcsimrootevent" )->SetAutoDelete(kTRUE);
+
+#ifdef mPMT
+ WCSimRootEvent * fHybridevent = new WCSimRootEvent();
+ // Set Branche
+ fInputTree->SetBranchAddress("wcsimrootevent2" ,&fHybridevent );
+ // Set autodelete to avoid memory leak
+ fInputTree->GetBranch("wcsimrootevent2" )->SetAutoDelete(kTRUE);
+#endif
+
+#ifdef OD_ON
+ WCSimRootEvent * fODevent = new WCSimRootEvent();
+ // Set Branche
+ fInputTree->SetBranchAddress("wcsimrootevent_OD",&fODevent );
+ // Set autodelete to avoid memory leak
+ fInputTree->GetBranch("wcsimrootevent_OD")->SetAutoDelete(kTRUE);
+#endif
+
+ // Read Geo
+ fInputGeoTree->GetEntry(0);
+
+
+ // Create Output TTree
+ TFile* fOutputFile = new TFile(sOutputFile.c_str(),"RECREATE");
+ fOutputFile->SetCompressionLevel(2);
+ TTree* fGeoTree = new TTree("wcsimGeoT","Geometry TTree");
+ TTree* fPrimaryTree = new TTree("Reduced","Reduced TTree");
+
+ // Set Branches
+ SetGeoBranch(fGeoTree);
+ SetCustomBranch(fPrimaryTree);
+
+ fGeoTree->Fill();
+
+ // Get PMT Number:
+ int nPMT_ID = fGeometry->GetWCNumPMT();
+#ifdef OD_ON
+ int nPMT_OD = fGeometry->GetODWCNumPMT();
+#endif
+ int nMultPMT = fGeometry->GetWCNumPMT(true);
+ // TODO solve this issue!
+
+ std::cout << " ID " << nPMT_ID << std::endl;
+ std::cout << " mPMT " << nMultPMT << std::endl;
+
+ HKManager::GetME()->SetGeometry(fGeometry,dDarkNoise * 1e3,dDarkNoiseHybrid * 1e3);
+
+ // Initialize LEAF
+ LEAF::GetME()->Initialize( HKManager::GetME()->GetGeometry(),
+ HKManager::GetME()->GetDarkNoise(),
+ HKManager::GetME()->GetGeometryPMT_ID(),
+ HKManager::GetME()->GetGeometryPMT_mPMT());
+ LEAF::GetME()->SetNThread(); // Set number of Threads, default in the class is 12
+
+ // Initialize HKAstroAnalysis
+#ifdef WITH_HK_ASTROANALYSIS
+ HKAstroAnalysis::GetME()->Initialize(HKManager::GetME()->GetGeometry());
+#endif
+
+#ifdef WITH_BONSAI
+ // Initialize Bonsai
+ WCSimBonsai* fBonsai = new WCSimBonsai();
+ fBonsai->Init(fGeometry);
+#endif
+ // Read Input Tree
+ int nPrimaryEvents = fInputTree->GetEntries();
+ int iWrite = 0;
+
+ TStopwatch timer;
+ timer.Reset();
+ timer.Start();
+
+
+ lf_hRelativeAngle = new TH1D("lf_ChargeProfile", "lf_ChargeProfile", 120, 0, 180);
+ lf_hRelativeAngleCos = new TH1D("lf_ChargeProfileCos", "lf_ChargeProfileCos", 100, -1, 1);
+ hRelativeAngle = new TH1D("ChargeProfile", "ChargeProfile", 120, 0, 180);
+ hRelativeAngleCos = new TH1D("ChargeProfileCos", "ChargeProfileCos", 100, -1, 1);
+ ToWall_Charge = new TH1D("ToWallCharge_hist", "ToWallCharge (Hist)", 1000, 0, 10000);
+ ToWall_Charge_2D = new TH2D("ToWallCharge", "ToWallCharge", 1000, 0, 10000, 50, 0, 400);
+ Sum_HitQ = new TH1D("Sum_HitQ", "Sum of HitQ per Wall bin", 1000, 0, 10000);
+ Count_Hits = new TH1D("Count_Hits", "Count of Hits per Wall bin", 1000, 0, 10000);
+
+ // Loop on Primary events
+ for(int i=StartEvent; i < StartEvent + NumberEvents; i++){
+ // Reset Hit vector
+ HKManager::GetME()->ResetHitInfo();
+
+ //if ( i%1000==0 ) {
+ timer.Stop();
+ std::cout << "Event # = " << i << " / " << nPrimaryEvents << " ( " << timer.RealTime() << " )\n";
+ timer.Reset();
+ timer.Start();
+ //}
+
+ fInputTree->GetEntry(i);
+
+ // Initialize output variables
+ eventId = iWrite;
+ triggerId = 0;
+ rawhit_num_noDN = 0;
+
+ true_particleId.clear();
+ true_energies.clear();
+ true_origins.clear();
+
+ digithit_pmtId.clear();
+ digithit_T.clear();
+ digithit_Q.clear();
+ Charge_PMT.clear();
+ digithit_dark.clear();
+ relativeAngle.clear();
+ lf_relativeAngle.clear();
+ rawhit_num = 0;
+ digithit_num = 0;
+
+ bs_vertex [0] = -9999.;
+ bs_vertex [1] = -9999.;
+ bs_vertex [2] = -9999.;
+ bs_vertex [3] = -9999.;
+ bs_good [0] = -9999.;
+ bs_good [1] = -9999.;
+ bs_good [2] = -9999.;
+
+ fBSTime = -9999;
+ fLFTime = -9999;
+
+ Hit_ID = 0;
+ Hit_ID_50 = 0;
+ Hit_ID_200 = 0;
+ Hit_ID_400 = 0;
+
+ Hit_mPMT = 0;
+ Hit_mPMT_50 = 0;
+ Hit_mPMT_200 = 0;
+ Hit_mPMT_400 = 0;
+
+ Hit_OD = 0;
+ Hit_OD_50 = 0;
+ Hit_OD_200 = 0;
+ Hit_OD_400 = 0;
+
+
+
+
+#ifdef WITH_BONSAI
+ // re-initialize Bonsai input
+ for ( int iHit = 0; iHit < 2000; iHit++ ) {
+ bsCAB[iHit] = 0;
+ bsT [iHit] = 0.;
+ bsQ [iHit] = 0.;
+ }
+ bsnhit[0] = 0;
+#endif
+ leaf_output.vtx = ROOT::Math::XYZTVector(0,0,0,0);
+ leaf_output.vtx_nll = -9999.;
+ //leaf_output.InTime = 0;
+ leaf_output.energy = 0.;
+ leaf_output.dir = ROOT::Math::XYZVector(0,0,0);
+ //leaf_output.SNRList = std::vector();
+
+ leaf_output_ana.Wall = 0;
+ for ( int iType = 0; iType < 3; iType++ ) {
+ leaf_output_ana.n50 [iType] = 0;
+ leaf_output_ana.dir [iType][0] = 0;
+ leaf_output_ana.dir [iType][1] = 0;
+ leaf_output_ana.dir [iType][2] = 0;
+ leaf_output_ana.dir_goodness [iType] = 0;
+ leaf_output_ana.dirKS [iType] = 0;
+ }
+
+ fLastRawHit = 0;
+
+ /****************************************************************************************/
+ /* ID events */
+ /****************************************************************************************/
+
+ //std::cout << " ID Event " << std::endl;
+ fHit = 0;
+ fHit_50 = 0;
+ fHit_200 = 0;
+ fHit_400 = 0;
+
+ /*bool bID =*/ AnalyseEvent(fIDevent,PMTType::kID,i);
+
+
+
+ Hit_ID = fHit;
+ Hit_ID_50 = fHit_50;
+ Hit_ID_200 = fHit_200;
+ Hit_ID_400 = fHit_400;
+
+ //std::cout << " Hit: " << Hit_ID << " event " << i << std::endl;
+
+ /****************************************************************************************/
+ /* mPMT events */
+ /****************************************************************************************/
+
+ //std::cout << " mPMT Event " << std::endl;
+ fHit = 0;
+ fHit_50 = 0;
+ fHit_200 = 0;
+ fHit_400 = 0;
+
+#ifdef mPMT
+ /*bool bmPMT =*/ AnalyseEvent(fHybridevent,PMTType::kmPMT,i);
+#endif
+ Hit_mPMT = fHit;
+ Hit_mPMT_50 = fHit_50;
+ Hit_mPMT_200 = fHit_200;
+ Hit_mPMT_400 = fHit_400;
+
+ //std::cout << " Hit: " << Hit_mPMT << " event " << i << std::endl;
+
+ /****************************************************************************************/
+ /* OD events */
+ /****************************************************************************************/
+#ifdef OD_ON
+ fHit = 0;
+ fHit_50 = 0;
+ fHit_200 = 0;
+ fHit_400 = 0;
+
+ // There should be a dedicated OD analyser, as many thing should be different than for ID
+ // Doesn't exist yet
+ /*bool bOD =*/ AnalyseODEvent(fODevent,PMTType::kOD);
+
+ Hit_OD = fHit;
+ Hit_OD_50 = fHit_50;
+ Hit_OD_200 = fHit_200;
+ Hit_OD_400 = fHit_400;
+#endif
+ /****************************************************************************************/
+ /* Benjamin Fitter */
+ /****************************************************************************************/
+
+ //std::cout << " Start LEAF " << std::endl;
+ TStopwatch timerLF;
+ timerLF.Reset();
+ timerLF.Start();
+
+ LEAF::GetME()->LoadHitsCollection(HKManager::GetME()->GetHitCollection_ID(),HKManager::GetME()->GetHitCollection_mPMT());
+ leaf_output = LEAF::GetME()->MakeSequentialFit();
+
+ timerLF.Stop();
+
+ fLFTime = timerLF.RealTime();
+
+ //std::cout << " n50 " << leaf_output_ana.n50[0] << " " << leaf_output_ana.n50[1] << " " << leaf_output_ana.n50[2] << std::endl;
+ std::cout << " LEAF took: " << timerLF.RealTime() << " for " << HKManager::GetME()->GetHitCollection_ID()->size() << " Hits"<< std::endl;
+ AnalyselfEvent(fIDevent,PMTType::kID);
+
+#ifdef WITH_HK_ASTROANALYSIS
+ HKAstroAnalysis::GetME()->SetVertex(leaf_output.Vtx);
+
+ leaf_output_ana.Wall = HKAstroAnalysis::GetME()->ComputeDistanceFromWall();
+
+ HKAstroAnalysis::GetME()->MakeAnalysis(HKManager::GetME()->GetHitCollection(),NormalPMT);
+
+ leaf_output_ana.n50 [NormalPMT] = HKAstroAnalysis::GetME()->Getn50();
+ leaf_output_ana.dirKS [NormalPMT] = HKAstroAnalysis::GetME()->GetdirKS();
+ leaf_output_ana.dir [NormalPMT][0] = HKAstroAnalysis::GetME()->Getdir_Simple()[0];
+ leaf_output_ana.dir [NormalPMT][1] = HKAstroAnalysis::GetME()->Getdir_Simple()[1];
+ leaf_output_ana.dir [NormalPMT][2] = HKAstroAnalysis::GetME()->Getdir_Simple()[2];
+ leaf_output_ana.dir_goodness [NormalPMT] = HKAstroAnalysis::GetME()->Getdir_Simple()[3];
+
+ HKAstroAnalysis::GetME()->MakeAnalysis(HKManager::GetME()->GetHitCollection(),MiniPMT);
+
+ leaf_output_ana.n50 [MiniPMT] = HKAstroAnalysis::GetME()->Getn50();
+ leaf_output_ana.dirKS [MiniPMT] = HKAstroAnalysis::GetME()->GetdirKS();
+ leaf_output_ana.dir [MiniPMT][0] = HKAstroAnalysis::GetME()->Getdir_Simple()[0];
+ leaf_output_ana.dir [MiniPMT][1] = HKAstroAnalysis::GetME()->Getdir_Simple()[1];
+ leaf_output_ana.dir [MiniPMT][2] = HKAstroAnalysis::GetME()->Getdir_Simple()[2];
+ leaf_output_ana.dir_goodness [MiniPMT] = HKAstroAnalysis::GetME()->Getdir_Simple()[3];
+
+ HKAstroAnalysis::GetME()->MakeAnalysis(HKManager::GetME()->GetHitCollection(),AllPMT);
+
+ leaf_output_ana.n50 [AllPMT] = HKAstroAnalysis::GetME()->Getn50();
+ leaf_output_ana.dirKS [AllPMT] = HKAstroAnalysis::GetME()->GetdirKS();
+ leaf_output_ana.dir [AllPMT][0] = HKAstroAnalysis::GetME()->Getdir_Simple()[0];
+ leaf_output_ana.dir [AllPMT][1] = HKAstroAnalysis::GetME()->Getdir_Simple()[1];
+ leaf_output_ana.dir [AllPMT][2] = HKAstroAnalysis::GetME()->Getdir_Simple()[2];
+ leaf_output_ana.dir_goodness [AllPMT] = HKAstroAnalysis::GetME()->Getdir_Simple()[3];
+
+ leaf_output_ana.Good = HKAstroAnalysis::GetME()->GoodnessBonsai();
+
+#endif
+
+ if (true_origins.size() > 0 ) {
+ ROOT::Math::XYZVector true_vertex(true_origins[0].X(), true_origins[0].Y(), true_origins[0].Z());
+ ROOT::Math::XYZVector reconstructed_vertex(leaf_output.vtx.X(), leaf_output.vtx.Y(), leaf_output.vtx.Z());
+
+ ROOT::Math::XYZVector lf_vertex(reconstructed_vertex);
+ ROOT::Math::XYZVector lf_Dire(leaf_output.dir);
+ lf_dWall = calculateDWall(R, h, lf_vertex);
+ lf_ToWall = calculateToWall(R, h, lf_vertex, lf_Dire);
+ lf_spatial_res = LEAFUtilities::GetDistance(true_vertex, reconstructed_vertex);
+ lf_jitter = abs(leaf_output.vtx.T()-true_origins[0].T());
+
+ TotalQ+=leaf_output.energy;
+ lf_Dir_interp = leaf_output.dir;
+ double angle_lf = lf_Dir_interp.Dot(particleDir);
+
+
+ lf_Dir_res = TMath::ACos(angle_lf)*180./TMath::Pi();
+
+
+ } else {
+ lf_spatial_res = -9999; // Assign invalid value if true vertex is missing
+ }
+#ifdef WITH_BONSAI
+ /****************************************************************************************/
+ /* Bonsai */
+ /****************************************************************************************/
+
+ if ( bsnhit[0] < 2000 && bsnhit[0] > 0 ) {
+
+ TStopwatch timerBS;
+ timerBS.Reset();
+ timerBS.Start();
+
+ // Some variable for Bonsai:
+ float bsvertex[4];
+ float bsresult[6];
+ float bsgood[3];
+ int bsnsel[1]; //nsel (SLE)
+
+ // Fit with Bonsai
+ //std::cout << "Bonsai hit: " << bsnhit[0] << std::endl;
+ //for ( int iHit = 0; iHit < bsnhit[0]; iHit++ ) {
+ // std::cout << " hit " << iHit << " cab " << bsCAB[iHit] << " T " << bsT[iHit] << " Q " << bsQ[iHit] << std::endl;
+ //}
+
+ fBonsai->BonsaiFit( bsvertex, bsresult, bsgood, bsnsel, bsnhit, bsCAB, bsT, bsQ);
+
+#ifdef WITH_HK_ASTROANALYSIS
+ HKAstroAnalysis::GetME()->SetVertex(bs_vertex);
+
+ bs_output_ana.Wall = HKAstroAnalysis::GetME()->ComputeDistanceFromWall();
+
+ HKAstroAnalysis::GetME()->MakeAnalysis(HKManager::GetME()->GetHitCollection());
+
+ bs_output_ana.n50 [NormalPMT] = HKAstroAnalysis::GetME()->Getn50();
+ bs_output_ana.dirKS [NormalPMT] = HKAstroAnalysis::GetME()->GetdirKS();
+ bs_output_ana.dir [NormalPMT][0] = HKAstroAnalysis::GetME()->Getdir_Simple()[0];
+ bs_output_ana.dir [NormalPMT][1] = HKAstroAnalysis::GetME()->Getdir_Simple()[1];
+ bs_output_ana.dir [NormalPMT][2] = HKAstroAnalysis::GetME()->Getdir_Simple()[2];
+ bs_output_ana.dir_goodness [NormalPMT] = HKAstroAnalysis::GetME()->Getdir_Simple()[3];
+
+ bs_output_ana.Good = HKAstroAnalysis::GetME()->GoodnessBonsai();
+#endif
+
+ bs_vertex[0] = bsvertex[0];
+ bs_vertex[1] = bsvertex[1];
+ bs_vertex[2] = bsvertex[2];
+ bs_vertex[3] = bsvertex[3];
+
+ bs_good[0] = bsgood[0];
+ bs_good[1] = bsgood[1];
+ bs_good[2] = bsgood[2];
+
+ //float diff = sqrt(pow(true_origin_X[0]- bsvertex[0],2.)+pow(true_origin_X[0]- bsvertex[0],2.)+pow(true_origin_X[0]- bsvertex[0],2.));
+
+ //timerBS.Stop();
+ //std::cout << " origin " << true_origin_X[0] << " " << true_origin_Y[0] << " " << true_origin_Z[0] << " diff: " << diff << std::endl;
+ //std::cout << " BS took: " << timerBS.RealTime() << " goodness " << bsgood[0]<< " " << bsgood[1]<< " " << bsgood[2] << " bsvertex: " << bsvertex[0] << " " << bsvertex[1] << " " << bsvertex[2] << std::endl;
+
+
+ fBSTime = timerBS.RealTime();
+ }
+#endif
+ /****************************************************************************************/
+ /* Fill output tree */
+ /****************************************************************************************/
+
+ fPrimaryTree->Fill();
+ iWrite += 1;
+ }
+ std::cout << "Event # = " << nPrimaryEvents << " / " << nPrimaryEvents << std::endl;
+// Convert histograms to graphs
+
+fOutputFile->cd();
+
+
+
+
+
+ lf_hRelativeAngle->Write();
+ lf_hRelativeAngleCos->Write();
+ hRelativeAngle->Write();
+ hRelativeAngleCos->Write();
+ ToWall_Charge->Write();
+ hToWall_RelAngle_Charge->Write();
+ pAngle_ToWall->Write();
+ hdWall_TotalCharge->Write();
+ hdWall_TotalCharge_corr->Write();
+ hdWall_TotalCharge_hist->Write();
+ ToWall_Charge_2D->Write();
+ hAngle_HitQ->Write();
+
+ fOutputFile->Write("",TObject::kOverwrite);
+
+ // Cleaning
+ delete fPrimaryTree;
+ delete fOutputFile;
+
+ //std::cout << "Total Charge: " << leaf_output.TotalCharge << std::endl;
+ return 1;
+//TotalQ 1000e10Mev: 207435 -- 100e10Mev: 20944.7 -- 1000e5MeV: 144119 -- 1000e3MeV: 45966.4
+}
+
+void SetCustomBranch(TTree* fPrimaryTree) {
+
+ fPrimaryTree->Branch("eventId", &eventId, "eventId/I");
+ fPrimaryTree->Branch("triggerId", &triggerId, "triggerId/I");
+
+ fPrimaryTree->Branch("true_particleId", &true_particleId);
+ fPrimaryTree->Branch("true_origins", &true_origins);
+ fPrimaryTree->Branch("true_energies", &true_energies);
+
+ fPrimaryTree->Branch("rawhit_num", &rawhit_num, "rawhit_num/I");
+ fPrimaryTree->Branch("digithit_num", &digithit_num, "digithit_num/I");
+
+ fPrimaryTree->Branch("ID_hits", &Hit_ID, "Hit_ID/I");
+ fPrimaryTree->Branch("ID_hits_50", &Hit_ID_50, "Hit_ID_50/I");
+ fPrimaryTree->Branch("ID_hits_200", &Hit_ID_200, "Hit_ID_200/I");
+ fPrimaryTree->Branch("ID_hits_400", &Hit_ID_400, "Hit_ID_400/I");
+
+ fPrimaryTree->Branch("mPMT_hits", &Hit_mPMT, "Hit_mPMT/I");
+ fPrimaryTree->Branch("mPMT_hits_50", &Hit_mPMT_50, "Hit_mPMT_50/I");
+ fPrimaryTree->Branch("mPMT_hits_200", &Hit_mPMT_200, "Hit_mPMT_200/I");
+ fPrimaryTree->Branch("mPMT_hits_400", &Hit_mPMT_400, "Hit_mPMT_400/I");
+
+ fPrimaryTree->Branch("bs_vertex", bs_vertex, "bs_vertex[4]/F");
+ fPrimaryTree->Branch("bs_good", bs_good, "bs_good[3]/F");
+ fPrimaryTree->Branch("bs_ctime", &fBSTime, "bs_ctime/F"); // Computation time
+
+ fPrimaryTree->Branch("lf_spatial_res", &lf_spatial_res, "lf_spatial_res/D");
+ fPrimaryTree->Branch("lf_jitter", &lf_jitter, "lf_jitter/D");
+ fPrimaryTree->Branch("DWall", &dWall, "DWall/D");
+ fPrimaryTree->Branch("lf_DWall", &lf_dWall, "DWall/D");
+ fPrimaryTree->Branch("lf_ToWall", &lf_ToWall, "DWall/D");
+
+ fPrimaryTree->Branch("reconstructed_energy", &leaf_output.energy, "reconstructed_energy/D");
+ fPrimaryTree->Branch("TotalCharge", &leaf_output.total_charge, "TotalCharge/D");
+ fPrimaryTree->Branch("ParticleDir", &particleDir,"particleDir[3]/D");
+ fPrimaryTree->Branch("lf_Dir", &leaf_output.dir);
+ fPrimaryTree->Branch("lf_Dir_res", &lf_Dir_res,"lf_Dir_res/D");
+ //fPrimaryTree->Branch("SNR", &leaf_output.SNRList);
+
+ fPrimaryTree->Branch("RelativeAngle", &relativeAngle);
+ fPrimaryTree->Branch("lf_RelativeAngle", &lf_relativeAngle);
+
+ fPrimaryTree->Branch("DigiHitQ", &digithit_Q);
+ fPrimaryTree->Branch("Charge_PMT", &Charge_PMT);
+ fPrimaryTree->Branch("DigiHitT", &digithit_T);
+
+
+
+
+ fPrimaryTree->Branch("lf_vertex", &leaf_output.vtx);
+ fPrimaryTree->Branch("lf_NLL", &leaf_output.vtx_nll, "lf_NLL/D");
+ //fPrimaryTree->Branch("lf_intime", &leaf_output.InTime, "lf_intime/I");
+ fPrimaryTree->Branch("lf_good", &leaf_output_ana.Good, "lf_good/D");
+ fPrimaryTree->Branch("lf_wall", &leaf_output_ana.Wall, "lf_wall/D");
+ fPrimaryTree->Branch("lf_n50", &leaf_output_ana.n50, "lf_n50[3]/I");
+ fPrimaryTree->Branch("lf_dir", &leaf_output_ana.dir, "lf_dir[3][3]/D");
+ fPrimaryTree->Branch("lf_dir_goodness", &leaf_output_ana.dir_goodness, "lf_dir_goodness[3]/D");
+ fPrimaryTree->Branch("lf_dirKS", &leaf_output_ana.dirKS, "lf_dirKS[3]/D");
+ fPrimaryTree->Branch("lf_ctime", &fLFTime, "lf_ctime/D"); // Computation time
+
+}
+
+void SetGeoBranch(TTree* fGeoTree) {
+
+ fGeoTree->Branch("wcsimrootgeom", fGeometry);
+
+}
+bool AnalyselfEvent(WCSimRootEvent * tEvent, PMTType pmtType){
+
+ int nVertex = 0; // Number of Vertex in event
+ int nTrack = 0; // Number of Track in event
+ int nRawCherenkovHits = 0; // Number of Raw Cherenkov hits
+ //int startingCherenkovHitID = 0; // starting ID of Digitized Cherenkov hits. Usually starts at 0, but if there are two types of PMTs, it does not.
+ int nDigitizedCherenkovHits = 0; // Number of Digitized Cherenkov hits
+ //double fTriggerTime = 0.;
+
+ // Declare some useful variables
+ WCSimRootTrigger * fRootTrigger;
+ //TClonesArray * fTimeArray;
+
+ // Currently only one Trigger is used (to be check)
+ //int iTrig = 0;
+ for(int iTrig = 0; iTrig < tEvent->GetNumberOfEvents(); iTrig++){
+
+ triggerId = iTrig;
+ fRootTrigger = tEvent->GetTrigger(iTrig);
+
+ // Grab the big arrays of times and parent IDs
+ //fTimeArray = fRootTrigger->GetCherenkovHitTimes();
+
+ // Get number of vertex and tracks
+#ifdef OLD_WCSIM
+ nVertex = 1;
+#else
+ nVertex = fRootTrigger->GetNvtxs();
+#endif
+
+ nTrack = fRootTrigger->GetNtrack();
+ nRawCherenkovHits = fRootTrigger->GetNumTubesHit();
+
+ if ( nTrack == 0 || nRawCherenkovHits == 0 ) {
+ // No track, no hit, nothing to do
+ return false;
+ }
+
+ // Get number of hits
+ //startingCherenkovHitID = digithit_pmtId.size();
+ nDigitizedCherenkovHits = fRootTrigger->GetNcherenkovdigihits();
+
+ if (nDigitizedCherenkovHits == 0){
+ continue;
+ }
+ //fTriggerTime = fRootTrigger->GetTriggerInfo()[2] - fRootTrigger->GetTriggerInfo()[1];
+ ROOT::Math::XYZVector Truepos;
+ // Loop on vertex
+ for(int iVertex = 0; iVertex < nVertex; iVertex++){
+#ifdef OLD_WCSIM
+ Truepos.SetXYZ(fRootTrigger->GetVtx(0),fRootTrigger->GetVtx(1),fRootTrigger->GetVtx(2));
+#else
+ Truepos.SetXYZ(fRootTrigger->GetVtxs(iVertex,0),fRootTrigger->GetVtxs(iVertex,1),fRootTrigger->GetVtxs(iVertex,2));
+#endif
+ if ( iVertex*2 > nTrack ) {
+ std::cout << "ERROR: Vertex and Track number incompatible (nVertex: " << nVertex << " ; nTrack " << nTrack << ") " << std::endl;
+ return 0;
+ }
+
+ // Beam info are registered in the Track
+ TObject * Track = (fRootTrigger->GetTracks())->At(iVertex*2);
+ WCSimRootTrack *wcTrack = dynamic_cast(Track);
+
+ if (wcTrack->GetParentId() == 0 && wcTrack->GetIpnu()==11) {
+ particleDir.SetXYZ(wcTrack->GetPdir(0),wcTrack->GetPdir(1),wcTrack->GetPdir(2));
+ particleDir = particleDir.Unit();
+ }
+ }
+
+ if ( nRawCherenkovHits < 1 ) return false;
+
+ std::vector times;
+
+ // Get number of hit
+ rawhit_num = nRawCherenkovHits;
+ // Loop on Digitized Hit
+ lf_Dir.SetXYZ(0,0,0);
+ double AllQ = 0;
+ for(int iDigitHit = 0; iDigitHit < nDigitizedCherenkovHits; iDigitHit++){
+ ROOT::Math::XYZVector lf_particleRelativePMTpos;
+ ROOT::Math::XYZVector particleRelativePMTpos;
+ ROOT::Math::XYZVector lf_relativePMTpos;
+ TObject *Hit = (fRootTrigger->GetCherenkovDigiHits())->At(iDigitHit);
+ WCSimRootCherenkovDigiHit *wcDigitHit = dynamic_cast(Hit);
+
+ int pmtId = wcDigitHit->GetTubeId();
+
+ //double HitT = wcDigitHit->GetT() + fTriggerTime;
+ double HitQ = wcDigitHit->GetQ();
+ int peForTube = wcDigitHit->GetQ();
+
+
+
+ WCSimRootPMT pmt;
+ pmt = fGeometry->GetPMT(pmtId - 1, false);
+ ROOT::Math::XYZVector PMTpos(pmt.GetPosition(0),pmt.GetPosition(1),pmt.GetPosition(2));
+ particleRelativePMTpos = (PMTpos - Truepos);
+ lf_relativePMTpos = (PMTpos - ROOT::Math::XYZVector(leaf_output.vtx.X(), leaf_output.vtx.Y(), leaf_output.vtx.Z()));
+ ROOT::Math::XYZVector vDir(particleRelativePMTpos.Unit());
+ ROOT::Math::XYZVector lfHitDir(lf_relativePMTpos.Unit());
+
+ lf_Dir += (lfHitDir * HitQ);
+
+ AllQ += HitQ;
+
+ double hitlfRelativeAngle = vDir.Dot(leaf_output.dir);
+ //std::cout << TMath::ACos(hitlfRelativeAngle)*180./TMath::Pi() << std::endl;
+
+ lf_relativeAngle.push_back((TMath::ACos(hitlfRelativeAngle))*180./TMath::Pi());
+
+ //std::cout << "Real Relative Angle : " << hitRelativeAngle*180./TMath::Pi() << std::endl;
+ //std::cout << "Charge of hit : " << HitQ << std::endl;
+ lf_hRelativeAngle->Fill((TMath::ACos(hitlfRelativeAngle))*180./TMath::Pi(), peForTube);
+ lf_hRelativeAngleCos->Fill(hitlfRelativeAngle, peForTube);
+ }
+ }
+ return true;
+
+}
+bool AnalyseEvent(WCSimRootEvent * tEvent, PMTType pmtType, int i) {
+
+ int nVertex = 0; // Number of Vertex in event
+ int nTrack = 0; // Number of Track in event
+ int nRawCherenkovHits = 0; // Number of Raw Cherenkov hits
+ int startingCherenkovHitID = 0; // starting ID of Digitized Cherenkov hits. Usually starts at 0, but if there are two types of PMTs, it does not.
+ int nDigitizedCherenkovHits = 0; // Number of Digitized Cherenkov hits
+ double fTriggerTime = 0.;
+
+ // Declare some useful variables
+ WCSimRootTrigger * fRootTrigger;
+ //TClonesArray * fTimeArray;
+
+ // Currently only one Trigger is used (to be check)
+ //int iTrig = 0;
+ for(int iTrig = 0; iTrig < tEvent->GetNumberOfEvents(); iTrig++){
+
+ triggerId = iTrig;
+ fRootTrigger = tEvent->GetTrigger(iTrig);
+
+ // Grab the big arrays of times and parent IDs
+ //fTimeArray = fRootTrigger->GetCherenkovHitTimes();
+
+ // Get number of vertex and tracks
+#ifdef OLD_WCSIM
+ nVertex = 1;
+#else
+ nVertex = fRootTrigger->GetNvtxs();
+#endif
+
+ nTrack = fRootTrigger->GetNtrack();
+ nRawCherenkovHits = fRootTrigger->GetNumTubesHit();
+
+ if ( nTrack == 0 || nRawCherenkovHits == 0 ) {
+ // No track, no hit, nothing to do
+ return false;
+ }
+
+ // Get number of hits
+
+ startingCherenkovHitID = digithit_pmtId.size();
+ nDigitizedCherenkovHits = fRootTrigger->GetNcherenkovdigihits();
+
+ if (nDigitizedCherenkovHits == 0){
+ continue;
+ }
+ fTriggerTime = fRootTrigger->GetTriggerInfo()[2] - fRootTrigger->GetTriggerInfo()[1];
+ ROOT::Math::XYZVector Truepos;
+ // Loop on vertex
+ for(int iVertex = 0; iVertex < nVertex; iVertex++){
+ // Beam info are registered in the Track
+ TObject * Track = (fRootTrigger->GetTracks())->At(iVertex*2);
+ WCSimRootTrack *wcTrack = dynamic_cast(Track);
+
+#ifdef OLD_WCSIM
+ Truepos.SetXYZ(fRootTrigger->GetVtx(0), fRootTrigger->GetVtx(1), fRootTrigger->GetVtx(2));
+ ROOT::Math::XYZTVector true_origin(fRootTrigger->GetVtx(0), fRootTrigger->GetVtx(1), fRootTrigger->GetVtx(2), wcTrack->GetTime());
+ true_origins.push_back(true_origin);
+
+#else
+ Truepos.SetXYZ(fRootTrigger->GetVtxs(iVertex,0), fRootTrigger->GetVtxs(iVertex,1), fRootTrigger->GetVtxs(iVertex,2));
+ ROOT::Math::XYZTVector true_origin(fRootTrigger->GetVtxs(iVertex,0), fRootTrigger->GetVtxs(iVertex,1), fRootTrigger->GetVtxs(iVertex,2), wcTrack->GetTime());
+ true_origins.push_back(true_origin);
+#endif
+ if ( iVertex*2 > nTrack ) {
+ std::cout << "ERROR: Vertex and Track number incompatible (nVertex: " << nVertex << " ; nTrack " << nTrack << ") " << std::endl;
+ return 0;
+ }
+
+
+ if (wcTrack->GetParentId() == 0 && wcTrack->GetIpnu()==11) {
+ particleDir.SetXYZ(wcTrack->GetPdir(0),wcTrack->GetPdir(1),wcTrack->GetPdir(2));
+ particleDir = particleDir.Unit();
+ }
+
+ leaf_output_ana.Wall = calculateToWall(R, h, Truepos, particleDir);
+ dWall = calculateDWall(R, h, Truepos);
+ std::cout << "true dir: " << particleDir.X() << "," << particleDir.Y() << "," << particleDir.Z() << std::endl;
+ std::cout << "true pos: " << Truepos.X() << ","<< Truepos.Y() << ","<< Truepos.Z() << std::endl;
+ //std::cout << "ToWall: " << leaf_output_ana.Wall << std::endl;
+ true_particleId.push_back( wcTrack->GetIpnu() );
+ true_energies.push_back( wcTrack->GetE() );
+
+ }
+
+ // Send True Position to LEAF for check
+ LEAF::GetME()->SetTrueVertexInfo(true_origins[0]);
+
+ if ( nRawCherenkovHits < 1 ) return false;
+
+ std::vector times;
+
+ // Get number of hit
+ rawhit_num = nRawCherenkovHits;
+ // Loop on Digitized Hit
+ lf_Dir.SetXYZ(0,0,0);
+ double AllQ = 0;
+
+
+
+ for(int iDigitHit = 0; iDigitHit < nDigitizedCherenkovHits; iDigitHit++){
+ ROOT::Math::XYZVector lf_particleRelativePMTpos;
+ ROOT::Math::XYZVector particleRelativePMTpos;
+ ROOT::Math::XYZVector lf_relativePMTpos;
+ TObject *Hit = (fRootTrigger->GetCherenkovDigiHits())->At(iDigitHit);
+ WCSimRootCherenkovDigiHit *wcDigitHit = dynamic_cast(Hit);
+
+ int pmtId = wcDigitHit->GetTubeId();
+
+ double HitT = wcDigitHit->GetT() + fTriggerTime;
+ double HitQ = wcDigitHit->GetQ();
+ int peForTube = wcDigitHit->GetQ();
+ AllQ += HitQ;
+ digithit_pmtId.push_back(pmtId);
+ digithit_T.push_back(HitT);
+ digithit_Q.push_back(HitQ);
+ Charge_PMT.push_back(peForTube);
+
+ Sum_HitQ->Fill(leaf_output_ana.Wall, HitQ);
+ Count_Hits->Fill(leaf_output_ana.Wall, 1);
+
+ Sum_HitQ_hdWall->Fill(dWall, HitQ);
+ Count_Hits_hdWall->Fill(dWall, 1);
+
+ WCSimRootPMT pmt;
+ pmt = fGeometry->GetPMT(pmtId - 1, false);
+
+ ROOT::Math::XYZVector PMTpos, PMTOrientation, vDir, lfHitDir;
+ double normPMTOrientation = 0, NormvDir = 0, NormlfHitDir = 0;
+
+ PMTpos.SetXYZ(pmt.GetPosition(0), pmt.GetPosition(1), pmt.GetPosition(2));
+
+ PMTOrientation.SetXYZ(pmt.GetOrientation(0), pmt.GetOrientation(1), pmt.GetOrientation(2));
+ PMTOrientation = PMTOrientation.Unit();
+
+ vDir = (PMTpos - Truepos).Unit();
+ lfHitDir = (PMTpos - ROOT::Math::XYZVector(leaf_output.vtx.X(),leaf_output.vtx.Y(),leaf_output.vtx.Z())).Unit();
+
+ normPMTOrientation = sqrt(normPMTOrientation);
+ NormvDir = sqrt(NormvDir);
+ NormlfHitDir = sqrt(NormlfHitDir);
+
+ double cosPMThitAngle = vDir.Dot(PMTOrientation);
+
+ double PMThitAngle = acos(-cosPMThitAngle) * (180.0 / M_PI);
+ //double cos2PMThitAngle = fabs(2 * cosPMThitAngle * cosPMThitAngle - 1);
+ //double Ftheta = 0.205 + 0.524*cos2PMThitAngle + 0.390*(pow(cos2PMThitAngle,2)) - 0.132*(pow(cos2PMThitAngle,3));
+
+ //std::cout << TMath::ACos(hitlfRelativeAngle)*180./TMath::Pi() << std::endl;
+ double hitRelativeAngle = vDir.Dot(particleDir);
+
+ relativeAngle.push_back((TMath::ACos(hitRelativeAngle))*180./TMath::Pi());
+ digithit_dark.push_back(false);
+ digithit_Type.push_back(pmtType);
+ //std::cout << "Real Relative Angle : " << hitRelativeAngle*180./TMath::Pi() << std::endl;
+ //std::cout << "Charge of hit : " << HitQ << std::endl;
+ hRelativeAngle->Fill((TMath::ACos(hitRelativeAngle))*180./TMath::Pi(), peForTube);
+ hRelativeAngleCos->Fill(hitRelativeAngle, peForTube);
+
+
+
+ /*WallPMTAngle_Charge->Fill(leaf_output_ana.Wall, PMThitAngle*2, HitQ)*/
+
+ hToWall_RelAngle_Charge->Fill(leaf_output_ana.Wall, (TMath::ACos(hitRelativeAngle))*180./TMath::Pi(), HitQ);
+ pAngle_ToWall->Fill((TMath::ACos(hitRelativeAngle))*180./TMath::Pi(),leaf_output_ana.Wall);
+ hAngle_HitQ->Fill(PMThitAngle, HitQ);
+
+ }
+
+
+ std::cout << "Hits: " << nDigitizedCherenkovHits << std::endl;
+ lf_Dir /= AllQ;
+ lf_Dir = lf_Dir.Unit();
+
+ for (int bin = 1; bin <= ToWall_Charge->GetNbinsX(); ++bin) {
+ double sum = Sum_HitQ->GetBinContent(bin);
+ double count = Count_Hits->GetBinContent(bin);
+
+ if (count > 0) {
+ ToWall_Charge->SetBinContent(bin, sum / count);
+ } else {
+ ToWall_Charge->SetBinContent(bin, 0);
+ }
+ }
+ ToWall_Charge_2D->Fill(leaf_output_ana.Wall, AllQ);
+
+ for (int bin = 1; bin <= hdWall_TotalCharge_hist->GetNbinsX(); ++bin) {
+ double sum = Sum_HitQ_hdWall->GetBinContent(bin);
+ double count = Count_Hits_hdWall->GetBinContent(bin);
+
+ if (count > 0) {
+ hdWall_TotalCharge_hist->SetBinContent(bin, sum / count);
+ } else {
+ hdWall_TotalCharge_hist->SetBinContent(bin, 0);
+ }
+ }
+
+ hdWall_TotalCharge->Fill(dWall, AllQ);
+ hdWall_TotalCharge_corr->Fill(dWall, leaf_output.energy);
+
+ /*double lf_rot_dir[3];
+ lf_rot_dir[0]=lf_Dir[2];
+ lf_rot_dir[1]=-lf_Dir[0];
+ lf_rot_dir[2]=-lf_Dir[1];
+ for(int j=0;j<3;j++){
+ lf_Dir[j] = lf_rot_dir[j];
+ }
+ lf_dir_res = sqrt(pow(lf_Dir[0] - particleDir[0], 2) + pow(lf_Dir[1] - particleDir[1], 2) + pow(lf_Dir[2] - particleDir[2], 2));*/
+
+ int iIdx_BS = 0;
+
+ nDigitizedCherenkovHits = digithit_pmtId.size();
+ digithit_num = nDigitizedCherenkovHits;
+
+ // Feed fitter
+ if ( pmtType == PMTType::kID ) {
+ for(int iDigitHit = startingCherenkovHitID; iDigitHit < nDigitizedCherenkovHits; iDigitHit++){
+ times.push_back(digithit_T[iDigitHit]);
+
+ //std::cout << " Add Hit with " << digithit_pmtId[iDigitHit] << " pmtType: " << pmtType << std::endl;
+ if ( digithit_pmtId[iDigitHit] <= 0 ) std::cout << " Weird PMT ID " << digithit_pmtId[iDigitHit] << std::endl;
+ HKManager::GetME()->AddHit_ID(digithit_T[iDigitHit],digithit_Q[iDigitHit],digithit_pmtId[iDigitHit]);
+
+ // Bonsai (do not store mPMT hits)
+ if ( iIdx_BS < 2000 && digithit_T[iDigitHit] < 4000. && digithit_T[iDigitHit] > -4000. ) {
+ //std::cout << pmtId << " " << HitT << std::endl;
+ bsCAB[iIdx_BS] = digithit_pmtId[iDigitHit];
+ bsT [iIdx_BS] = digithit_T[iDigitHit]; // shift BS time is needed if interaction time is < 0, this needs to be considered
+ bsQ [iIdx_BS] = digithit_Q[iDigitHit];
+
+ iIdx_BS += 1;
+ }
+ }
+ }
+ else {
+ for(int iDigitHit = startingCherenkovHitID; iDigitHit < nDigitizedCherenkovHits; iDigitHit++){
+ times.push_back(digithit_T[iDigitHit]);
+ //std::cout << " Add Hit with " << digithit_pmtId[iDigitHit] << " pmtType: " << pmtType << std::endl;
+ if ( digithit_pmtId[iDigitHit] <= 0 ) std::cout << " Weird PMT ID " << digithit_pmtId[iDigitHit] << std::endl;
+ HKManager::GetME()->AddHit_mPMT(digithit_T[iDigitHit],digithit_Q[iDigitHit],digithit_pmtId[iDigitHit]);
+ }
+ }
+
+ // Bonsai
+ bsnhit[0] = iIdx_BS;
+ fHit = digithit_pmtId.size();
+
+ // Compute hits
+ std::vector Hit_time_50;
+ std::vector Hit_time_200;
+ std::vector Hit_time_400;
+
+ std::sort(times.begin(),times.end());
+ for(int iDigitHit = 0; iDigitHit < nDigitizedCherenkovHits; iDigitHit++){
+
+ double HitT = times[iDigitHit];
+
+ Hit_time_50.push_back(HitT);
+ Hit_time_200.push_back(HitT);
+ Hit_time_400.push_back(HitT);
+
+ // Count hit in 50 ns window
+ while ( HitT - Hit_time_50[0] > 50. )
+ Hit_time_50.erase(Hit_time_50.begin());
+
+ // Count hit in 200 ns window
+ while ( HitT - Hit_time_200[0] > 200. )
+ Hit_time_200.erase(Hit_time_200.begin());
+
+
+ // Count hit in 400 ns window
+ while ( HitT - Hit_time_400[0] > 400. )
+ Hit_time_400.erase(Hit_time_400.begin());
+
+ if ( (unsigned int) fHit_50 < Hit_time_50.size() ) fHit_50 = Hit_time_50.size();
+ if ( (unsigned int) fHit_200 < Hit_time_200.size() ) fHit_200 = Hit_time_200.size();
+ if ( (unsigned int) fHit_400 < Hit_time_400.size() ) fHit_400 = Hit_time_400.size();
+ }
+
+ }
+
+
+ return true;
+}
\ No newline at end of file
diff --git a/example/env.sh b/example/env.sh
new file mode 100644
index 0000000..d388c56
--- /dev/null
+++ b/example/env.sh
@@ -0,0 +1,6 @@
+export LEAFDIR=../
+export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$LEAFDIR/lib
+export ROOT_INCLUDE_PATH=${LEAFDIR}/leaf/DataModel-lite
+
+export WCSIMDIR=${LEAFDIR}/libWCSIM
+export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$WCSIMDIR
\ No newline at end of file
diff --git a/example/test_example.sh b/example/test_example.sh
index 9a21b8b..63c37b2 100755
--- a/example/test_example.sh
+++ b/example/test_example.sh
@@ -1,16 +1,38 @@
#!/bin/bash
+# Launch single job
+#
+#SBATCH --job-name=leaf # Job name
+#SBATCH --output=/sps/t2k/lperisse/Soft/leaf/logs/analysis_%j.log # Standard output and error log
+#SBATCH --partition=htc # Partition choice
+#SBATCH --ntasks=1 # Maximum number of parallel processes
+#SBATCH --mem=5G # Amount of memory required
+#SBATCH --time=0-01:00 # 7 days by default on htc partition
+
+
# Replace $input contents with the path toward your WCSim input file
# (Note: make sure the WCSim path in RunAtStart matches the WCSim version you used to produce the file)
-input=$WCSIMDIR/wcsim_hybrid.root
+input=/sps/t2k/lperisse/Soft/wcsim/results/electron/wcsim1p12p29_electron_HK_test.root
+output=/sps/t2k/lperisse/Soft/wcsim/results/electron/wcsim1p12p29_electron_HK_test.leaf.root
+
-# test.root doesn't exist, this is just a placeholder
+# Argument syntaxe to run ./analysis:
+#-i Input WCSim ROOT file
+#-o Output ROOT file
+# d Dark noise frequency of PMTs in kHz
+# h Dark noise frequency of mPMTS in kHz
+# s First event to analyze
+#-e Last event to analyze
+# v VERBOSE option
if [ ! -f "$input" ]; then
echo "Set input file doesn't exist! Add your WCSim input file here. (Note: make sure the WCSim path in RunAtStart matches the WCSim version you used to produce the file)"
else
echo "Processing ${input}..."
- ./analysis $input output.root 0 0
+ cd ${LEAFDIR}/example
+ ./analysis -i $input -o $output -s 0 -e 5 -d 4.2 -h 0.0 -v
fi
+
+
diff --git a/inputs/.gitignore b/inputs/.gitignore
new file mode 100644
index 0000000..e69de29
diff --git a/inputs/PDF_electron_upto100MeV_uniform_isotropic_10k_withDN.root b/inputs/PDF_electron_upto100MeV_uniform_isotropic_10k_withDN.root
new file mode 100644
index 0000000..c834078
Binary files /dev/null and b/inputs/PDF_electron_upto100MeV_uniform_isotropic_10k_withDN.root differ
diff --git a/inputs/timePDF_DRnew_Large.root b/inputs/timePDF_DRnew_Large.root
deleted file mode 100644
index 0b1b78c..0000000
Binary files a/inputs/timePDF_DRnew_Large.root and /dev/null differ
diff --git a/inputs/timePDF_Directionality.root b/inputs/timePDF_Directionality.root
deleted file mode 100644
index 2baa242..0000000
Binary files a/inputs/timePDF_Directionality.root and /dev/null differ
diff --git a/inputs/timePDF_Directionality_DRnew.root b/inputs/timePDF_Directionality_DRnew.root
deleted file mode 100644
index 0289350..0000000
Binary files a/inputs/timePDF_Directionality_DRnew.root and /dev/null differ
diff --git a/inputs/timePDF_HE.root b/inputs/timePDF_HE.root
deleted file mode 100644
index f81ca2c..0000000
Binary files a/inputs/timePDF_HE.root and /dev/null differ
diff --git a/leaf/DataModel-lite/Constants/Mass.hpp b/leaf/DataModel-lite/Constants/Mass.hpp
new file mode 100644
index 0000000..35a1878
--- /dev/null
+++ b/leaf/DataModel-lite/Constants/Mass.hpp
@@ -0,0 +1,13 @@
+#ifndef MASS_HPP
+#define MASS_HPP
+
+#include