diff --git a/.gitignore b/.gitignore index 4f2b247..c4d3cce 100644 --- a/.gitignore +++ b/.gitignore @@ -1,8 +1,12 @@ *.root *.o *.swp +*.root +*.log .depend leaf/DataModelRootDict* leaf/HKAstroAnalysis.* leaf/DataModel lib/* +job.sh +Config.sh \ No newline at end of file diff --git a/example/ConfigExample.sh b/example/ConfigExample.sh new file mode 100644 index 0000000..273dfbf --- /dev/null +++ b/example/ConfigExample.sh @@ -0,0 +1,32 @@ +#!/bin/bash +x + +# This script is used to set up the environment for running LEAF + +export HitTimeLimitsNegative=-4 +export HitTimeLimitsPositive=8 +export SearchVtxStep=300 +export SearchVtxTolerance=20 +export maxHitsAngle=190 +export N_Neighbors=5 +export maxDistanceToNeighbors=6800 + +export STimePDFLimitsQueueNegative=-3; +export STimePDFLimitsQueuePositive=4; +export STimePDFLimitsQueueNegative_fullTimeWindow=0; +export STimePDFLimitsQueuePositive_fullTimeWindow=0; +export TimeWindowSizeFull=1500; +export Averaging=20; # Number of points we used to average the vertex on. + +export MinimizeLimitsNegative=-700; +export MinimizeLimitsPositive=1000; + +export IntegrationTimeWindow=50; + +export DirTolerance=30; + +export StepByStep=false; # Step by step mode was not test and is not supported by multithreading +export UseDirectionality=false; +export HighEnergy=false; +export DoubleNLL=true; +export Limit_mPMT=true; +export DirTakeAll=false; \ No newline at end of file diff --git a/example/GNUmakefile b/example/GNUmakefile index 3d2baf2..ede518e 100644 --- a/example/GNUmakefile +++ b/example/GNUmakefile @@ -10,19 +10,23 @@ 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$(BONSAIDIR)/bonsai +INCFLAGS += -I$(WCSIM_BUILD_DIR)/include/WCSim +# INCFLAGS += -I$(WCSIMDIR)/include +# INCFLAGS += -I$(BONSAIDIR)/bonsai INCFLAGS += -I$(LEAFDIR)/leaf INCFLAGS += -I$(LEAFDIR)/leaf/DataModel -LIBS += -L${WCSIMDIR} -lWCSimRoot +# LIBS += -L${WCSIMDIR} -lWCSimRoot +LIBS += -L${WCSIM_BUILD_DIR}/lib -lWCSimRoot LIBS += $(shell root-config --libs) -lMinuit -LIBS_BONSAI += -L${BONSAIDIR} -lWCSimBonsai +# LIBS_BONSAI += -L${BONSAIDIR} -lWCSimBonsai LIBS_LEAF += -L${LEAFDIR}/lib -lDataModelLite -lHKManager -lLEAF #-lHKAstroAnalysis OBJECT = analysis +CXXFLAGS += -std=c++17 + all: $(OBJECT) analysis: analysis.o diff --git a/example/analysis.cpp b/example/analysis.cpp index 0534ec8..2f62396 100644 --- a/example/analysis.cpp +++ b/example/analysis.cpp @@ -4,775 +4,1176 @@ /** Date: Febuary 10th 2020 **/ /** Desc: Example application code for Benjamin's Low-E Fitter for Hyper-K **/ /*********************************************************************************/ +#include "analysis.h" -#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.hh" -#include "HKManager.hh" -//#define WITH_BONSAI -#ifdef WITH_BONSAI -#include "WCSimBonsai.hh" -#endif +///* UTILITIES -#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 -// HK_ASTROANALYSIS is an analysis class computing some informations like n50, or the bonsai-like goodness -// Ask G. Pronost to have access to this class (you need to be a member of SK or HK collaboration) -//#define WITH_HK_ASTROANALYSIS +// Apply analysis on the event -#ifdef WITH_HK_ASTROANALYSIS -#include "HKAstroAnalysis.hh" -#endif +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; -#define ID_EVENT 1 -#define OD_EVENT 2 -#define mPMT_EVENT 3 -#define UNDEFINED_EVENT 0 + 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; - int eventId; // event Id - int triggerId; // trigger Id (with event Id) + 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; - - - struct FitterAnalysis { - double Wall; - double Good; - int n50[3]; - double dir[3][3]; - double dir_goodness[3]; - - double dirKS[3]; - }; - + return t_wall; +} - //True information - 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) - - //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 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 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}); +} - LEAF::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; - - WCSimRootGeom * fGeometry; - - // Input variable for Bonsai - int bsCAB[2000]; - float bsT[2000]; - float bsQ[2000]; - int bsnhit[1]; +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."); + } -void SetCustomBranch(TTree* fPrimaryTree); -void SetCustomBranchInput(TTree* fPrimaryTree); -void SetGeoBranch(TTree* fGeoTree); + // Cosine of the angle + double cosTheta = dot(a,b) / (magA * magB); -// Apply analysis on the event -bool AnalyseEvent(WCSimRootEvent * tEvent, int iEventType); + // 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; +} -int main(int argc, char** argv){ +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; +} + +// void InitializeHistograms() +// { + // 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); + + + + // Hit_Q->GetXaxis()->SetTitle("Hit Charge [p.e.]"); + // Hit_Q->GetYaxis()->SetTitle("Number of Hits"); +// } + + +///* MAIN EXECUTION + +int main(int argc, char **argv) +{ + // InitializeHistograms(); std::string sInputFile = ""; std::string sOutputFile = ""; - int iNeededArgc = 3; - double dDarkNoise = 0.; //kHz - double dDarkNoiseHybrid = 0.; //kHz + failAmount = 0; + + int iNeededArgc = 3; + double dDarkNoise = 0.; // kHz + double dDarkNoiseHybrid = 0.; // kHz + iNeededArgc += 1; -#ifdef mPMT +#ifdef mPMT iNeededArgc += 1; #endif - - if ( argc == iNeededArgc ) { + + 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; + 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 { - + else + { + std::cout << "Synthax: " << argv[0] << " input output"; std::cout << " DN_in_kHz_B&L"; - -#ifdef mPMT + +#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"); - + 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(); + TTree *fInputTree = (TTree *)fInputFile->Get("wcsimT"); + TTree *fInputGeoTree = (TTree *)fInputFile->Get("wcsimGeoT"); + + fLeafGeometry = 0; + fInputGeoTree->SetBranchAddress("wcsimrootgeom", &fLeafGeometry); + + WCSimRootEvent *fIDevent = new WCSimRootEvent(); // Set Branche - fInputTree->SetBranchAddress("wcsimrootevent" ,&fIDevent ); + fInputTree->SetBranchAddress("wcsimrootevent", &fIDevent); // Set autodelete to avoid memory leak - fInputTree->GetBranch("wcsimrootevent" )->SetAutoDelete(kTRUE); - - WCSimRootEvent * fHybridevent = new WCSimRootEvent(); + fInputTree->GetBranch("wcsimrootevent")->SetAutoDelete(kTRUE); + + // WCSimRootEvent *fHybridevent = new WCSimRootEvent(); #ifdef mPMT // Set Branche - fInputTree->SetBranchAddress("wcsimrootevent2" ,&fHybridevent ); + // fInputTree->SetBranchAddress("wcsimrootevent2", &fHybridevent); // Set autodelete to avoid memory leak - fInputTree->GetBranch("wcsimrootevent2" )->SetAutoDelete(kTRUE); + fInputTree->GetBranch("wcsimrootevent2")->SetAutoDelete(kTRUE); #endif - + #ifdef OD_ON - WCSimRootEvent * fODevent = new WCSimRootEvent(); + WCSimRootEvent *fODevent = new WCSimRootEvent(); // Set Branche - fInputTree->SetBranchAddress("wcsimrootevent_OD",&fODevent ); + 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"); + 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"); - + TTree *fGeoTree = new TTree("wcsimGeoT", "Geometry TTree"); + TTree *fPrimaryTree = new TTree("Reduced", "Reduced TTree"); + // Set Branches SetGeoBranch(fGeoTree); - SetCustomBranch(fPrimaryTree); - + + FitterOutput leaf_output; + SetCustomBranch(fPrimaryTree, leaf_output); + fGeoTree->Fill(); - + // Get PMT Number: - int nPMT_ID = fGeometry->GetWCNumPMT(); + int nPMT_ID = fLeafGeometry->GetWCNumPMT(); #ifdef OD_ON - int nPMT_OD = fGeometry->GetODWCNumPMT(); + int nPMT_OD = fLeafGeometry->GetODWCNumPMT(); #endif - int nMultPMT = fGeometry->GetWCNumPMT(true); - // TODO solve this issue! - + int nMultPMT = fLeafGeometry->GetWCNumPMT(true); + std::cout << " ID " << nPMT_ID << std::endl; std::cout << " mPMT " << nMultPMT << std::endl; - - HKManager::GetME()->SetGeometry(fGeometry,dDarkNoise * 1e3,dDarkNoiseHybrid * 1e3); - + + HKManager::GetME()->SetGeometry(fLeafGeometry, dDarkNoise * 1e3, dDarkNoiseHybrid * 1e3); + // Initialize LEAF LEAF::GetME()->Initialize(HKManager::GetME()->GetGeometry()); - 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 + + 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) + { + 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; + + envValue = std::getenv("N_Neighbors"); + 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; + + envValue = std::getenv("maxDistanceToNeighbors"); + 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; + // Read Input Tree int nPrimaryEvents = fInputTree->GetEntries(); - int iWrite = 0; - + + if(maxEvents > 0) nPrimaryEvents = std::min(nPrimaryEvents, maxEvents); + + int iWrite = 0; + TStopwatch timer; timer.Reset(); timer.Start(); - - // Loop on Primary events - for(int i=0; i < nPrimaryEvents; i++){ - + + // Loop on Primary events + for (int i = 0; i < nPrimaryEvents; 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(); + // HKManager::GetME()->ResetSecondaryHitInfo(); + + // if ( i%1000==0 ) { + timer.Stop(); + std::cout << "Event # = " << i << " / " << nPrimaryEvents << " ( " << timer.RealTime() << " )\n"; + timer.Reset(); + timer.Start(); //} - - fInputTree->GetEntry(i); - + + fInputTree->GetEntry(i); + // Initialize output variables - eventId = iWrite; - triggerId = 0; - rawhit_num_noDN = 0; - + eventId = iWrite; + nTrigger = 0; + usedTriggerId = 0; + 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(); + 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_dark.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; + digithit_Angle.clear(); + digithit_NormAngle.clear(); + digithit_NeighborsDist.clear(); + hit_is_DR.clear(); + Charge_PMT.clear(); + relativeAngle.clear(); + lf_relativeAngle.clear(); + hit_residual.clear(); + // digithit_Angle_NLL.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_5_15ns = 0; + hit_50ns = 0; + + Hit_ID = 0; + Hit_ID_50 = 0; + Hit_ID_200 = 0; + Hit_ID_400 = 0; + + Hit_mPMT = 0; + Hit_mPMT_20 = 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[0] = 0.; - leaf_output.Vtx[1] = 0.; - leaf_output.Vtx[2] = 0.; - leaf_output.Vtx[3] = 0.; - leaf_output.NLL = -9999.; - leaf_output.InTime = 0; - - 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; + + Hit_OD = 0; + Hit_OD_50 = 0; + Hit_OD_200 = 0; + Hit_OD_400 = 0; + + // leaf_output.Vtx[0] = 0.; + // leaf_output.Vtx[1] = 0.; + // leaf_output.Vtx[2] = 0.; + // leaf_output.Vtx[3] = 0.; + // leaf_output.NLL = -9999.; + // leaf_output.InTime = 0; + // leaf_output.Energy = 0.; + // leaf_output.Dir[0] = 0.; + // leaf_output.Dir[1] = 0.; + // leaf_output.Dir[2] = 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; - + + truepos = std::vector(4); + trueDir = std::vector(3); + lf_Dir = std::vector(3); + + // for (int j = 0; j < 3; j++) + // { + // trueDir[j] = 0.; + // lf_Dir[j] = 0.; + // lf_Dir_interp[j] = 0.; + // } + lf_spatial_res = 0.; + lf_Dir_res = 0.; + lf_time_res = 0.; + lf_energy_res = 0.; + + bestTrigger = 0; + fLfTriggerTime = 0.; + /****************************************************************************************/ /* ID events */ /****************************************************************************************/ - - //std::cout << " ID Event " << std::endl; - fHit = 0; - fHit_50 = 0; + + // std::cout << " ID Event " << std::endl; + fHit = 0; + fHit_20 = 0; + fHit_50 = 0; fHit_200 = 0; fHit_400 = 0; - /*bool bID =*/ AnalyseEvent(fIDevent,ID_EVENT); + bool shouldProcess = AnalyseEvent(fIDevent, ID_EVENT); - Hit_ID = fHit; - Hit_ID_50 = fHit_50; + if(!shouldProcess) + { + failAmount++; + std::cout << "Event # = " << i << " cannot be processed " << std::endl; + continue; + } + + Hit_ID = fHit; + Hit_ID_20 = fHit_20; + Hit_ID_50 = fHit_50; Hit_ID_200 = fHit_200; Hit_ID_400 = fHit_400; - //std::cout << " Hit: " << Hit_ID << " event " << i << std::endl; - + // std::cout << " Hit: " << Hit_ID << " event " << i << std::endl; + /****************************************************************************************/ /* mPMT events */ /****************************************************************************************/ - - //std::cout << " mPMT Event " << std::endl; - fHit = 0; - fHit_50 = 0; + + // std::cout << " mPMT Event " << std::endl; + fHit = 0; + fHit_20 = 0; + fHit_50 = 0; fHit_200 = 0; fHit_400 = 0; - + #ifdef mPMT - /*bool bmPMT =*/ AnalyseEvent(fHybridevent,mPMT_EVENT); + // /*bool bmPMT =*/AnalyseEvent(fHybridevent, mPMT_EVENT); #endif - Hit_mPMT = fHit; - Hit_mPMT_50 = fHit_50; + Hit_mPMT = fHit; + Hit_mPMT_20 = fHit_20; + Hit_mPMT_50 = fHit_50; Hit_mPMT_200 = fHit_200; Hit_mPMT_400 = fHit_400; - - //std::cout << " Hit: " << Hit_mPMT << " event " << i << std::endl; - + + // std::cout << " Hit: " << Hit_mPMT << " event " << i << std::endl; + /****************************************************************************************/ /* OD events */ /****************************************************************************************/ + #ifdef OD_ON - fHit = 0; - fHit_50 = 0; + fHit = 0; + fHit_20 = 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,OD_EVENT); + /*bool bOD =*/AnalyseODEvent(fODevent, OD_EVENT); - Hit_OD = fHit; - Hit_OD_50 = fHit_50; + Hit_OD = fHit; + Hit_OD_50 = fHit_50; Hit_OD_200 = fHit_200; Hit_OD_400 = fHit_400; -#endif +#endif /****************************************************************************************/ /* Benjamin Fitter */ /****************************************************************************************/ - //std::cout << " Start LEAF " << std::endl; + // std::cout << " Start LEAF " << std::endl; TStopwatch timerLF; timerLF.Reset(); timerLF.Start(); - + // To be replaced by meaningful TriggerTime TimeDelta fDummyTrigger(0.); - leaf_output = LEAF::GetME()->MakeFit(HKManager::GetME()->GetHitCollection(),fDummyTrigger); - -#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 - + + leaf_output = LEAF::GetME()->MakeSequentialFit(HKManager::GetME()->GetHitCollection(), fDummyTrigger); + // leaf_output = LEAF::GetME()->MakeJointFit(HKManager::GetME()->GetHitCollection(), fDummyTrigger); + + // stepOneHasTrueVtx = leaf_output.stepOneContainsTrueVtx; + // Vtx_Search_ComputeTime = leaf_output.Vtx_Search_ComputeTime; + // Vtx_Minimize_ComputeTime = leaf_output.Vtx_Minimize_ComputeTime; + // goodness = leaf_output.NLLR; + 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()->Size() << " Hits"<< std::endl; -#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 + // 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()->Size() << " Hits"; + std::cout << " (Vtx Search : " << leaf_output.Vtx_Search_ComputeTime << " , Vtx Minimize : " << leaf_output.Vtx_Minimize_ComputeTime << ", Dir Search : " << leaf_output.Dir_Search_ComputeTime << " , Dir Minimize : " << leaf_output.Dir_Minimize_ComputeTime << " , Energy Fit : " << leaf_output.Energy_Fit_ComputeTime << ")" << std::endl; + bool validEvent = PostLeafAnalysis(fIDevent,ID_EVENT, leaf_output); + + // std::cout << " After LEAF " << std::endl; /****************************************************************************************/ /* Fill output tree */ /****************************************************************************************/ - - fPrimaryTree->Fill(); - iWrite += 1; + if(validEvent) + { + std::cout << " output vertex at last : " << leaf_output.Vtx[0] << " " << leaf_output.Vtx[1] << " " << leaf_output.Vtx[2] << " " << leaf_output.Vtx[3] << std::endl; + fPrimaryTree->Fill(); + iWrite += 1; + } } - + std::cout << "Event # = " << nPrimaryEvents << " / " << nPrimaryEvents << std::endl; - fOutputFile->Write("",TObject::kOverwrite); - - // Cleaning + std::cout << ((nPrimaryEvents - failAmount) / nPrimaryEvents) * 100 << "% of the events were processed" << std::endl; + + fOutputFile->cd(); + + fOutputFile->Write("", TObject::kOverwrite); + delete fPrimaryTree; delete fOutputFile; - return 1; + // gInterpreter->GenerateDictionary("vector", "vector"); + + // TFile *fReInputFile = TFile::Open(sOutputFile.c_str()); + // if (!fReInputFile || fReInputFile->IsZombie()) { + // std::cerr << "Error: Could not open output file " << sOutputFile << std::endl; + // } else { + // TTree *reTree = nullptr; + // fReInputFile->GetObject("Reduced", reTree); + // if (!reTree) { + // std::cerr << "Error: Could not find TTree 'Reduced' in file " << sOutputFile << std::endl; + // } else { + // std::vector *ReVertex = nullptr; + // reTree->SetBranchAddress("lf_vertex", &ReVertex); + // if (reTree->GetEntries() > 0) { + // reTree->GetEntry(0); + // if (ReVertex && ReVertex->size() >= 4) { + // std::cout << "reVertex : " << ReVertex->at(0) << " " << ReVertex->at(1) << " " << ReVertex->at(2) << " " << ReVertex->at(3) << std::endl; + // } else { + // std::cerr << "Error: 'lf_vertex' branch is missing or does not have enough elements." << std::endl; + // } + // } else { + // std::cerr << "Error: TTree 'Reduced' has no entries." << std::endl; + // } + // } + // fReInputFile->Close(); + // delete fReInputFile; + // } + + return 1; } -void SetCustomBranch(TTree* fPrimaryTree) { +void SetCustomBranch(TTree *fPrimaryTree, FitterOutput leaf_output) +{ - fPrimaryTree->Branch("eventId", &eventId, "eventId/I"); - fPrimaryTree->Branch("triggerId", &triggerId, "triggerId/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("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_vertex", &leaf_output.Vtx, "lf_vertex[4]/D"); - fPrimaryTree->Branch("lf_NLL", &leaf_output.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 + fPrimaryTree->Branch("eventId", &eventId, "eventId/I"); + fPrimaryTree->Branch("nTrigger", &nTrigger, "nTrigger/I"); + fPrimaryTree->Branch("triggerUsed", &usedTriggerId, "triggerUsed/I"); + fPrimaryTree->Branch("rawhit_num", &rawhit_num, "rawhit_num/I"); + 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("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_20", &Hit_mPMT_20, "Hit_mPMT_20/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("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_intime", &leaf_output.InTime, "lf_intime/I"); + fPrimaryTree->Branch("lf_good", &leaf_output.NLLR, "lf_good/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 + 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_MyDir_DNLL", &leaf_output.myDNLL, "lf_MyDir_DNLL/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"); -void SetGeoBranch(TTree* fGeoTree) { + fPrimaryTree->Branch("stepOneHasTrueVtx", &leaf_output.stepOneContainsTrueVtx, "stepOneHasTrueVtx/I"); + fPrimaryTree->Branch("stepOneHasTrueDir", &leaf_output.stepOneContainsTrueDir, "stepOneHasTrueDir/I"); - fGeoTree->Branch("wcsimrootgeom", fGeometry); + fPrimaryTree->Branch("Leaf_ComputeTime", &leaf_output.Leaf_ComputeTime, "Leaf_ComputeTime/D"); + fPrimaryTree->Branch("Vtx_Search_ComputeTime", &leaf_output.Vtx_Search_ComputeTime, "Vtx_Search_ComputeTime/D"); + fPrimaryTree->Branch("Vtx_Minimize_ComputeTime", &leaf_output.Vtx_Minimize_ComputeTime, "Vtx_Minimize_ComputeTime/D"); + fPrimaryTree->Branch("Dir_Quick_Search_ComputeTime", &leaf_output.Dir_Quick_Search_ComputeTime, "Dir_Quick_Search_ComputeTime/D"); + fPrimaryTree->Branch("Dir_Search_ComputeTime", &leaf_output.Dir_Search_ComputeTime, "Dir_Search_ComputeTime/D"); + fPrimaryTree->Branch("Dir_Minimize_ComputeTime", &leaf_output.Dir_Minimize_ComputeTime, "Dir_Minimize_ComputeTime/D"); + fPrimaryTree->Branch("Energy_Fit_ComputeTime", &leaf_output.Energy_Fit_ComputeTime, "Energy_Fit_ComputeTime/D"); + + fPrimaryTree->Branch("rawTriggerTime", &rawTriggerTime, "rawTriggerTime/D"); + fPrimaryTree->Branch("TotalCharge", &leaf_output.TotalCharge, "TotalCharge/D"); + fPrimaryTree->Branch("SNR", &leaf_output.SNRList); + fPrimaryTree->Branch("RelativeAngle", &relativeAngle); + fPrimaryTree->Branch("lf_RelativeAngle", &lf_relativeAngle); + 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); + + // fPrimaryTree->Branch("digithit_Angle_NLL", &digithit_Angle_NLL); + fPrimaryTree->Branch("EnergyRatio", &EnergyRatio,"EnergyRatio/D"); } -bool AnalyseEvent(WCSimRootEvent * tEvent, int iEventType) { +void SetGeoBranch(TTree *fGeoTree) +{ + fGeoTree->Branch("wcsimrootgeom", fLeafGeometry); +} +bool AnalyseEvent(WCSimRootEvent *tEvent, int iEventType) +{ int iHybrid = 0; - if ( iEventType == mPMT_EVENT ) iHybrid = 1; + // if ( iEventType == mPMT_EVENT ) iHybrid = 1; + 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 nDigitizedCherenkovHits = 0; // Number of Digitized Cherenkov hits + nTrigger = tEvent->GetNumberOfEvents(); - 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++){ + usedTriggerId = -1; + rawTriggerTime = 0.; + WCSimRootTrigger *fRootTrigger; + + TimeDelta fDummyTrigger(0.); + TimeDelta fTimeCorrection = HKManager::GetME()->GetHitCollection()->timestamp - fDummyTrigger; //*still have no idea what this is - triggerId = iTrig; + //* find best trigger + int bestNumbOfHits = 0; + bestTrigger = -1; + for (int iTrig = 0; iTrig < nTrigger; iTrig++) + { fRootTrigger = tEvent->GetTrigger(iTrig); - - // Grab the big arrays of times and parent IDs - //fTimeArray = fRootTrigger->GetCherenkovHitTimes(); - - // Get number of vertex and tracks + nDigitizedCherenkovHits = fRootTrigger->GetNcherenkovdigihits(); + + if(nDigitizedCherenkovHits > bestNumbOfHits && fRootTrigger->GetTriggerInfo().size() >= 3) + { + bestNumbOfHits = nDigitizedCherenkovHits; + bestTrigger = iTrig; + } + } + + if(bestTrigger < 0) return false; //* Cannot process this event, none of the triggers have sufficent information + + usedTriggerId = bestTrigger; + + //* Retreive true infos (onlyh set by wcsim in the first trigger) + + for (int iTrig = 0; iTrig < tEvent->GetNumberOfEvents(); iTrig++) + { + // std::cout << "number of hits : " << tEvent->GetTrigger(iTrig)->GetNcherenkovhits() << std::endl; + // nTrigger = iTrig; + fRootTrigger = tEvent->GetTrigger(iTrig); + #ifdef OLD_WCSIM - nVertex = 1; + nVertex = 1; #else - nVertex = fRootTrigger->GetNvtxs(); + 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(); - fTriggerTime = fRootTrigger->GetTriggerInfo()[2] - fRootTrigger->GetTriggerInfo()[1]; - - // Loop on vertex - for(int iVertex = 0; iVertex < nVertex; iVertex++){ - + nTrack = fRootTrigger->GetNtrack(); + nRawCherenkovHits = fRootTrigger->GetNumTubesHit(); + + nDigitizedCherenkovHits = fRootTrigger->GetNcherenkovdigihits(); + + //* Fist trigger is the one used to compute true Vertex & Dir + if (iTrig == 0) + { + // Loop on vertex + for (int iVertex = 0; iVertex < nVertex; iVertex++) + { + #ifdef OLD_WCSIM - true_origin_X.push_back( fRootTrigger->GetVtx(0) ); - true_origin_Y.push_back( fRootTrigger->GetVtx(1) ); - true_origin_Z.push_back( fRootTrigger->GetVtx(2) ); + for(int j=0; j<3; j++) truepos[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)); #else - 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) ); -#endif - if ( iVertex*2 > nTrack ) { - std::cout << "ERROR: Vertex and Track number incompatible (nVertex: " << nVertex << " ; nTrack " << nTrack << ") " << std::endl; - return 0; + for(int j=0; j<3; j++) truepos[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)); +#endif + + // LEAF::GetME()->SetTrueVertexInfo(true_origin_X[0],true_origin_Y[0],true_origin_Z[0],0); + if (iVertex * 2 > nTrack) + { + std::cout << "ERROR: Vertex and Track number incompatible (nVertex: " << nVertex << " ; nTrack " << nTrack << ") " << std::endl; + continue; + } + + // 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) + { + for(int j=0; j<3; j++) trueDir[j] = wcTrack->GetPdir(j); + Normalize(trueDir); + // double magnitude=sqrt(pow(trueDir[0], 2) + pow(trueDir[1], 2) + pow(trueDir[2], 2)); + // for(int j=0; j<3; j++) trueDir[j] = trueDir[j]/magnitude; + } + + toWall = calculateToWall(R, h, truepos, trueDir); + dWall = calculateDWall(R, h, truepos); + // std::cout << "true dir: " << trueDir[0] << "," << trueDir[1] << "," << trueDir[2] << std::endl; + // std::cout << "true pos: " << truepos[0] << ","<< truepos[1] << ","<< truepos[2] << std::endl; + + true_particleId.push_back(wcTrack->GetIpnu()); + true_energy.push_back(wcTrack->GetE()); + true_origin_T.push_back(fRootTrigger->GetVtx(3)); + // true_origin_T.push_back(wcTrack->GetTime()); + + // std::cout << "true origin T = " << true_origin_T[0] << std::endl; + // std::cout << "true origin xyz " << true_origin_X[0] << " " << true_origin_Y[0] << " " << true_origin_Z[0] << std::endl; } - - // Beam info are registered in the Track - TObject * Track = (fRootTrigger->GetTracks())->At(iVertex*2); - WCSimRootTrack *wcTrack = dynamic_cast(Track); - - true_particleId.push_back( wcTrack->GetIpnu() ); - true_energy.push_back( wcTrack->GetE() ); - true_origin_T.push_back( wcTrack->GetTime() ); + + // Send True Position to LEAF for check + std::vector vVtxTrue(3, 0.); + vVtxTrue[0] = true_origin_X[0]; + vVtxTrue[1] = true_origin_Y[0]; + vVtxTrue[2] = true_origin_Z[0]; + // vVtxTrue[3] = true_origin_T[0]; + + // std::vector trueVtx = {true_origin_X[0], true_origin_Y[0], true_origin_Z[0]}; + // std::cout << "True Vertex: " << vVtxTrue[0] << " " << vVtxTrue[1] << " " << vVtxTrue[2] << std::endl; + SetTrueVertexInfo(vVtxTrue, true_origin_T[0]); + SetTrueDirInfo(trueDir); + // std::cout << "Set True Vertex done" << std::endl; + + rawhit_num = nRawCherenkovHits; } - - // Send True Position to LEAF for check - std::vector vVtxTrue(3,0.); - vVtxTrue[0] = true_origin_X[0]; - vVtxTrue[1] = true_origin_Y[0]; - vVtxTrue[2] = true_origin_Z[0]; - - //LEAF::GetME()->SetTrueVertexInfo(vVtxTrue,0); - - if ( nRawCherenkovHits < 1 ) return false; - + + if(iTrig != usedTriggerId) continue; //* The trigger being processed + + auto triggerInfo = fRootTrigger->GetTriggerInfo(); + fLfTriggerTime = triggerInfo[2] - triggerInfo[1]; // time of th event - offset (common offset of all events, depends on wcsim config giles) + rawTriggerTime = triggerInfo[2]; + + startingCherenkovHitID = digithit_pmtId.size(); std::vector times; - - // Get number of hit - rawhit_num = nRawCherenkovHits; + + //* 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++){ + for (int iDigitHit = 0; iDigitHit < nDigitizedCherenkovHits; iDigitHit++) + { + //* Get and Compute Hit Informations TObject *Hit = (fRootTrigger->GetCherenkovDigiHits())->At(iDigitHit); - WCSimRootCherenkovDigiHit *wcDigitHit = dynamic_cast(Hit); - - int pmtId = wcDigitHit->GetTubeId(); + 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(); + int peForTube = wcDigitHit->GetQ(); + // double PMTpos[3]; + std::vector lfHitDir = std::vector(3); + std::vector PMTOrientation = std::vector(3); + std::vector vDir = std::vector(3); + // double normPMTOrientation = 0; + // double NormvDir = 0, NormlfHitDir = 0; + bool isDR = IsDarkRateHit(fRootTrigger, wcDigitHit); + double hitResidual = GetResidualTime(truepos, true_origin_T[0], wcDigitHit, fTimeCorrection, fLfTriggerTime); + for (int j = 0; j < 3; j++) + { + // PMTpos[j] = pmt.GetPosition(j); + PMTOrientation[j] = pmt.GetOrientation(j); + // normPMTOrientation += PMTOrientation[j] * PMTOrientation[j]; + vDir[j] = pmt.GetPosition(j) - truepos[j]; + // lfHitDir[j] = PMTpos[j] - leaf_output.Vtx[j]; + lfHitDir[j] = pmt.GetPosition(j); + // NormvDir += vDir[j] * vDir[j]; + // NormlfHitDir += lfHitDir[j] * lfHitDir[j]; + } + // normPMTOrientation = sqrt(normPMTOrientation); + // NormvDir = sqrt(NormvDir); + // NormlfHitDir = sqrt(NormlfHitDir); - double HitT = wcDigitHit->GetT() + fTriggerTime; - double HitQ = wcDigitHit->GetQ(); + Normalize(vDir); + Normalize(lfHitDir); + Normalize(PMTOrientation); + + std::vector neg_vDir = { -vDir[0], -vDir[1], -vDir[2] }; + double NormAngle = dot(neg_vDir, PMTOrientation); + + // double cosPMThitAngle = 0; + for (int j = 0; j < 3; j++) lf_Dir[j] += lfHitDir[j] * HitQ; + // { + // PMTOrientation[j] /= normPMTOrientation; + // vDir[j] /= NormvDir; + // lfHitDir[j] /= NormlfHitDir; + // cosPMThitAngle += vDir[j] * PMTOrientation[j]; + // } + // 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[0]*trueDir[0]+vDir[1]*trueDir[1]+vDir[2]*trueDir[2]); + double hitRelativeAngle = dot(vDir, trueDir); + + // 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 distance = Astro_GetDistance(PMTpos, truepos); + // double tof = distance / fLightSpeed; + // double hitTime = (fTimeCorrection + HitT) / TimeDelta::ns; + // double hitResidual = hitTime - tof - true_origin_T[0]; + + //* 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; + + //* Filter base on distance to hit neighbors + double distanceToNeighbor = GetDistanceToNeighbors(fRootTrigger, wcDigitHit, N_Neighbors); + // if(distanceToNeighbor > maxDistanceToNeighbors) continue; + + //*DR FILTER + // if(isDR) continue; //! HERE + + // std::cout<< "IS DR : " << isDR << std::endl; + + ///* fill infos to output root file + + //* Tree Branches + 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); - digithit_dark.push_back(false); + 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); } - + + for(long unsigned int j=0;j -5) hit_5_15ns++; + if(hit_residual[j] < 50 && hit_residual[j] > -50) hit_50ns++; + } + + // std::cout << hit_50ns << " hits in 50 ns window, " << hit_5_15ns << " hits in 5-15 ns window "; + + for(int j=0;j<3;j++) lf_Dir[j] /= AllQ; + Normalize(lf_Dir); + + // double NormlfDir = TMath::Sqrt(lf_Dir[0]*lf_Dir[0]+lf_Dir[1]*lf_Dir[1]+lf_Dir[2]*lf_Dir[2]); + + // for(int j=0;j<3;j++) lf_Dir[j] /= NormlfDir; + int iIdx_BS = 0; nDigitizedCherenkovHits = digithit_pmtId.size(); - digithit_num = nDigitizedCherenkovHits; + digithit_num = nDigitizedCherenkovHits - startingCherenkovHitID; + + // std::cout << digithit_num << " digitized hits in total" << std::endl; // Feed fitter - for(int iDigitHit = startingCherenkovHitID; iDigitHit < nDigitizedCherenkovHits; iDigitHit++){ - - //iDigitHit+=startingCherenkovHitID; + for (int iDigitHit = startingCherenkovHitID; iDigitHit < nDigitizedCherenkovHits; iDigitHit++) + { + + // iDigitHit+=startingCherenkovHitID; times.push_back(digithit_T[iDigitHit]); - //std::cout << " Add Hit with " << digithit_pmtId[iDigitHit] << " Hybrid: " << iHybrid << std::endl; - 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]); + // std::cout << " Add Hit with " << digithit_pmtId[iDigitHit] << " Hybrid: " << iHybrid << std::endl; + 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]); + // if(iTrig == bestTrigger) + // { + // HKManager::GetME()->AddSecondaryHit(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; + 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]; - + 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; } } - + // Bonsai - bsnhit[0] = iIdx_BS; - fHit = digithit_pmtId.size(); - + bsnhit[0] = iIdx_BS; + fHit = digithit_pmtId.size(); + // Compute hits + std::vector Hit_time_20; 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++){ - + + std::sort(times.begin(), times.end()); + for (int iDigitHit = 0; iDigitHit < nDigitizedCherenkovHits; iDigitHit++) + { double HitT = times[iDigitHit]; - + + Hit_time_20.push_back(HitT); Hit_time_50.push_back(HitT); Hit_time_200.push_back(HitT); Hit_time_400.push_back(HitT); - + + // Count hit in 20 ns window + while (HitT - Hit_time_20[0] > 20.) Hit_time_20.erase(Hit_time_20.begin()); + // 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()); - + while (HitT - Hit_time_50[0] > 50.) Hit_time_50.erase(Hit_time_50.begin()); + + // Count hit in 200 ns windowAnalyseEvent + 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(); + while (HitT - Hit_time_400[0] > 400.) Hit_time_400.erase(Hit_time_400.begin()); + + if ((unsigned int)fHit_20 < Hit_time_20.size()) fHit_20 = Hit_time_20.size(); + 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; +} + +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]); + + std::cout << " output vertex : " << leaf_output.Vtx[0] << " " << leaf_output.Vtx[1] << " " << leaf_output.Vtx[2] << " " << leaf_output.Vtx[3] << std::endl; + + // double correctionSum = 0; + // for (long unsigned int i = 0; i < leaf_output.AngleCorrections.size(); i++) + // { + // corrections += leaf_output.AngleCorrections[i]; + // correctionAmount+=1; + // } + + totalChargeMean += leaf_output.Energy / leaf_output.TotalCharge; + totalChargeAmounts += 1; + EnergyRatio = totalChargeMean / totalChargeAmounts; + + // std::cout << " Mean correction : " << corrections / correctionAmount << std::endl; + // std::cout << " Total charge : " << leaf_output.TotalCharge << " and energy : " << leaf_output.Energy << std::endl; + // std::cout << " Mean charge ratio " << totalChargeMean / totalChargeAmounts << std::endl; + + + //* RESOLUTIONS COMPUTATION + std::vector true_vertex = {true_origin_X[0], true_origin_Y[0], true_origin_Z[0]}; + // std::vector reconstructed_vertex = {leaf_output.Vtx[0], leaf_output.Vtx[1], leaf_output.Vtx[2]}; + + // double lf_vertex[3]={reconstructed_vertex[0], reconstructed_vertex[1], reconstructed_vertex[2]}; + // double lf_Dire[3]={leaf_output.Dir[0],leaf_output.Dir[1],leaf_output.Dir[2]}; + 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]); + + // if(!(lf_dWall < 600 && lf_ToWall < 850 && lf_ToWall <= 40 + lf_dWall * 2)) return false; + + // digithit_Angle_NLL = leaf_output.VtxHitHangles; - } + // TotalQ+=leaf_output.Energy; + // double lf_Dir_interp[3]; + // lf_Dir_interp[0] = leaf_output.Dir[0]; + // lf_Dir_interp[1] = leaf_output.Dir[1]; + // lf_Dir_interp[2] = leaf_output.Dir[2]; + // double angle_lf = lf_Dir_interp[0]*trueDir[0]+lf_Dir_interp[1]*trueDir[1]+lf_Dir_interp[2]*trueDir[2]; - return true; -} + // double angle_lf = dot(leaf_output.Dir, trueDir); + + //* ANGLE ANALYSIS + WCSimRootTrigger * fRootTrigger = tEvent->GetTrigger(bestTrigger); + // TimeDelta fDummyTrigger(0.); + // TimeDelta fTimeCorrection = HKManager::GetME()->GetHitCollection()->timestamp - fDummyTrigger; + // int nDigitizedCherenkovHits = fRootTrigger->GetNcherenkovdigihits(); + // if (nDigitizedCherenkovHits == 0) continue; + + 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); + // int pmtId = wcDigitHit->GetTubeId(); + // double HitQ = wcDigitHit->GetQ(); + // int peForTube = wcDigitHit->GetQ(); + // WCSimRootPMT pmt; + // double PMTpos[3]; + + // bool isDR = IsDarkRateHit(fRootTrigger, wcDigitHit); + + // double hitResidual = GetResidualTime(truepos, true_origin_T[0], wcDigitHit, fTimeCorrection, fLfTriggerTime); + + 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) - truepos[j]; + Normalize(trueHitDir); + // for(int j=0;j<3;j++) + // { + // PMTpos[j] = pmt.GetPosition(j); + //PMTOrientation[j] = pmt.Orientation(j); + // trueHitDir[j] = pmt.GetPosition(j) - truepos[j]; + // lfHitDir[j] = pmt.GetPosition(j) - leaf_output.Vtx[j]; + // } + + // std::vector vDir = std::vector(3); + // double lfHitDir[3]; + + // for(int j=0;j<3;j++) + // { + // vDir[j] = particleRelativePMTpos[j]; + // lfHitDir[j] = lf_relativePMTpos[j]; + // } + //Normalize(PMTOrientation); + + for(int j=0;j<3;j++) lf_Dir[j] += lfHitDir[j]*wcDigitHit->GetQ(); + + // double hitlfRelativeAngle = (vDir[0]*leaf_output.Dir[0]+vDir[1]*leaf_output.Dir[1]+vDir[2]*leaf_output.Dir[2]); + double hitlfRelativeAngle = dot(trueHitDir, leaf_output.Dir); + //std::cout << TMath::ACos(hitlfRelativeAngle)*180./TMath::Pi() << std::endl; + + //std::cout << "Real Relative Angle : " << hitRelativeAngle*180./TMath::Pi() << std::endl; + //std::cout << "Charge of hit : " << HitQ << std::endl; + 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.h new file mode 100644 index 0000000..df73f78 --- /dev/null +++ b/example/analysis.h @@ -0,0 +1,203 @@ + +#include +#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.hh" +#include "LeafDefinitions.hh" +#include "HKManager.hh" + +#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 + +//-----------------------------------------------------------------------------------------// + +int eventId; // event Id +int nTrigger; // trigger Id (with event Id) +int usedTriggerId; // trigger Id used for the fit + +// double TotalQ; + +// double trueDir[3] = {0, 0, 0}; + +struct FitterAnalysis +{ + // double Wall; + // double Good; + int n50[3]; + double dir[3][3]; + double dir_goodness[3]; + + double dirKS[3]; +}; + +// True information +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 leaf_Vertex; // the estimated vertex +std::vector leaf_Dir; // the estimated direction +std::vector leaf_MyDir; // direction with one prior that is the the quick direction +std::vector leaf_QuickDir; // quick computed directions form to hit vectors + +std::vector relativeAngle; +std::vector lf_relativeAngle; +std::vector lf_Dir; +std::vector trueDir; +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 correctedDigithit_T; +std::vector hit_is_DR; +std::vector digithit_Q; // List of hit charge' +std::vector digithit_Angle; // List of hit charge' +std::vector digithit_Angle_NLL; // todo delete this +std::vector digithit_NormAngle; // List of hit charge' +std::vector digithit_NeighborsDist; // List of hit charge' +std::vector Charge_PMT; +std::vector hit_residual; // List of residual times +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_time_res; +double lf_Dir_res; +double lf_Quick_Dir_res; +double lf_MyDir_res; +double lf_energy_res; + +// FitterOutput leaf_output; +FitterAnalysis leaf_output_ana; +FitterAnalysis bs_output_ana; + +int Hit_ID; +int Hit_ID_20; +int Hit_ID_50; +int Hit_ID_200; +int Hit_ID_400; + +int Hit_mPMT; +int Hit_mPMT_20; +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 hit_5_15ns; //between -5 and 15ns +int hit_50ns; + +int fHit = 0; +int fHit_20 = 0; +int fHit_50 = 0; +int fHit_200 = 0; +int fHit_400 = 0; + +double dWall =0; +double toWall =0; +double lf_dWall =0; +double lf_ToWall =0; +// double goodness; + +double R = 3242.96; +double h = 6701.41; + +double maxHitAngle = 190.0; +double N_Neighbors = 5; +double maxDistanceToNeighbors = 8000.0; + +double corrections; +double correctionAmount; +double totalChargeMean; +double totalChargeAmounts; +double EnergyRatio; + +WCSimRootGeom *fLeafGeometry; + +// Input variable for Bonsai +int bsCAB[2000]; +float bsT[2000]; +float bsQ[2000]; +int bsnhit[1]; + +int failAmount; +int startingCherenkovHitID = 0; // starting ID of Digitized Cherenkov hits. Usually starts at 0, but if there are two types of PMTs, it does not. +std::vector truepos; + +double fLfTriggerTime; +int bestTrigger; + +bool failed; + +// int stepOneHasTrueVtx; +// double Vtx_Search_ComputeTime; +// double Vtx_Minimize_ComputeTime; +double rawTriggerTime; + +//-----------------------------------------------------------------------------------------// + +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); diff --git a/leaf/DataModel-lite/CreateDataModel-lite.sh b/leaf/DataModel-lite/CreateDataModel-lite.sh index 3482c06..e645a0b 100755 --- a/leaf/DataModel-lite/CreateDataModel-lite.sh +++ b/leaf/DataModel-lite/CreateDataModel-lite.sh @@ -12,7 +12,6 @@ else rsync -av --exclude=*.o $HK_ASTROANALYSIS_DIR/DataModelRoot/PMTInfo.* . rsync -av --exclude=*.o $HK_ASTROANALYSIS_DIR/DataModelRoot/TimeDelta.* . rsync -av --exclude=*.o $HK_ASTROANALYSIS_DIR/DataModelRoot/Pos3DT.* . - echo "You should update manually Environments.h" fi diff --git a/leaf/GNUmakefile b/leaf/GNUmakefile index be73889..e8194e0 100755 --- a/leaf/GNUmakefile +++ b/leaf/GNUmakefile @@ -8,7 +8,6 @@ HEADERS = LEAF.hh include ../Makefile/Makefile.${OSNAME} -ROOT_VERSION = $(shell root-config --version | tr -d '/' ) # set compiler options for ROOT @@ -18,20 +17,40 @@ CXXFLAGS_BASE += '-fPIC' -Wall -Wpedantic -Wno-long-long -DDATAMODEL_LITE CXXFLAGS = $(CXXFLAGS_BASE) #CXXFLAGS_ROOT = $(CXXFLAGS_BASE) -ifeq ($(shell echo "$(ROOT_VERSION) >= 6.0" | bc -l ), 1) +# #ifeq ($(shell echo "$(ROOT_VERSION) >= 6.0" | bc -l ), 1) +# ROOT_VERSION = $(shell root-config --version | tr -d '/' ) +# ifeq ($(shell [ $(ROOT_VERSION) -ge 6.0 ] && echo 1 || echo 0), 1) +# MAKE_LEAF_WITH_ROOT6 = 1 +# CXXFLAGS += -std=c++17 -DROOT6 # modifies to c++17 (Nicolas) +# #CXXFLAGS_ROOT += -std=c++11 -DROOT6 +# else +# CXXFLAGS += -std=c++11 -DROOT5 +# endif + +# ROOT_VERSION = $(shell root-config --version | tr -d '/') +# ifeq ($(shell echo "$(ROOT_VERSION) >= 6.0" | bc -l), 1) +# MAKE_LEAF_WITH_ROOT6 = 1 +# CXXFLAGS += -std=c++17 -DROOT6 # modifies to c++17 (Nicolas) +# #CXXFLAGS_ROOT += -std=c++11 -DROOT6 +# else +# CXXFLAGS += -std=c++11 -DROOT5 +# endif + + +ROOT_VERSION = $(shell root-config --version | cut -d'.' -f1,2) +ifeq ($(shell echo "$(ROOT_VERSION) >= 6.0" | bc -l), 1) MAKE_LEAF_WITH_ROOT6 = 1 -CXXFLAGS += -std=c++11 -DROOT6 +CXXFLAGS += -std=c++17 -DROOT6 # modifies to c++17 (Nicolas) #CXXFLAGS_ROOT += -std=c++11 -DROOT6 else -CXXFLAGS += -std=c++11 -DROOT5 -#CXXFLAGS_ROOT += -DROOT5 +CXXFLAGS += -std=c++17 -DROOT5 endif INCFLAGS = -I. -I$(shell root-config --incdir) -INCFLAGS += -I$(WCSIMDIR)/include +INCFLAGS += -I$(WCSIM_BUILD_DIR)/include/WCSim INCFLAGS += -I./DataModel -LIBS += -L${WCSIMDIR} -lWCSimRoot +LIBS += -L${WCSIM_BUILD_DIR}/lib -lWCSimRoot LIBS += $(shell root-config --libs) -lMinuit LIB_DATAMODEL = ../lib/libDataModelLite.so @@ -60,7 +79,7 @@ else endif # library -$(LIB_FITTER): LEAF.o +$(LIB_FITTER): LeafConfig.o LeafInputs.o LeafUtility.o LeafSplines.o Likelihoods.o LEAF.o @if [ ! -d ../lib ]; then \ mkdir ../lib; \ fi diff --git a/leaf/HKManager.hh b/leaf/HKManager.hh index 7584da8..4703f54 100644 --- a/leaf/HKManager.hh +++ b/leaf/HKManager.hh @@ -33,7 +33,7 @@ //#define WCSIM_single_PMT_type // Number of PMT configuration: -#define NPMT_CONFIGURATION 2 +#define NPMT_CONFIGURATION 1 class WCSimReader { @@ -45,9 +45,11 @@ class WCSimReader { void SetGeometry( WCSimRootGeom * wGeo, double dDarkRate_Normal=4200., double dDarkRate_mPMT=100. ); const HitCollection* GetHitCollection() { return &fHitCollection; } + // const HitCollection* GetSecondaryHitCollection() { return &fSecondaryHitCollection; } // Manage hit info void ResetHitInfo() { fHitCollection.Clean(); } + // void ResetSecondaryHitInfo() { fSecondaryHitCollection.Clean(); } void AddHit(double time, double charge, int pmtType, int tubeNumber) { // tubeNumber is from 1 to xxx in Hit array @@ -59,6 +61,18 @@ class WCSimReader { fHitCollection.Add(hHit); } + + // void AddSecondaryHit(double time, double charge, int pmtType, int tubeNumber) { + + // // tubeNumber is from 1 to xxx in Hit array + // if ( tubeNumber < 1 ) { + // std::cout << "ERROR: tubeNumber is below 0 (" << tubeNumber << ")" << std::endl; + // } + + // Hit hHit (tubeNumber, time, charge, (HKAA::PMTType) pmtType); + + // fSecondaryHitCollection.Add(hHit); + // } /* // Commented on 2020/11/19 by Guillaume: I don't remember their purpose. To be removed? @@ -76,6 +90,7 @@ class WCSimReader { // HitCollection HitCollection fHitCollection; + HitCollection fSecondaryHitCollection; // Darknoise double fDarkRate_Normal; diff --git a/leaf/LEAF.cc b/leaf/LEAF.cc index c2aca75..d2ade06 100644 --- a/leaf/LEAF.cc +++ b/leaf/LEAF.cc @@ -9,976 +9,512 @@ #include "LEAF.hh" #include "TStopwatch.h" -double LEAF::fSTimePDFLimitsQueueNegative = 0; -double LEAF::fSTimePDFLimitsQueuePositive = 0; -LEAF* LEAF::myFitter=NULL; - +LEAF *LEAF::myFitter = NULL; std::mutex mtx; /*****************************************************************************************************/ -void MinuitLikelihood(int& /*nDim*/, double * /*gout*/, double & NLL, double par[], int /*flg*/){ - //TStopwatch timer; - //timer.Reset(); - //timer.Start(); - - std::vector vertexPosition(4,0.);//In centimeters - for(int i=0;i<4;i++) vertexPosition[i]=par[i]; - //int pmtType=par[4]; - int nhits=par[5]; - double lowerLimit=par[6]; - double upperLimit=par[7]; - //double expoSigma=par[8]; - double directionality=par[9]; - - NLL=LEAF::GetME()->FindNLL_Likelihood(vertexPosition,nhits,lowerLimit,upperLimit,true,false,directionality); - - //timer.Stop();7 - //std::cout << " Likelihood: " << timer.RealTime() << std::endl; - -} - +// void MinimizeVertex_CallThread(int iStart, int iIte, std::vector> initialVertex, double *limits, double stepSize, int nhits, int nCandidates, int tolerance, int verbose, bool likelihood, bool average, double lowerLimit, double upperLimit, int directionality) +// { +// LEAF::GetME()->MinimizeVertex_thread(iStart, iIte, initialVertex, limits, stepSize, nhits, nCandidates, tolerance, verbose, likelihood, average, lowerLimit, upperLimit, directionality); +// } -void MinimizeVertex_CallThread( - int iStart, int iIte, - std::vector< std::vector > initialVertex, double * limits, double stepSize, int nhits, - int nCandidates, int tolerance, int verbose,bool likelihood,bool average,double lowerLimit, double upperLimit, int directionality){ - - LEAF::GetME()->MinimizeVertex_thread(iStart, iIte, - initialVertex, limits, stepSize, nhits, - nCandidates, tolerance, verbose, likelihood, average, lowerLimit, upperLimit, directionality); - -} - -void SearchVertex_CallThread( - int iStart, int iIte, - int nhits,int tolerance,bool likelihood,double lowerLimit, double upperLimit,int directionality){ - - LEAF::GetME()->SearchVertex_thread(iStart, iIte, - nhits, tolerance, likelihood, lowerLimit, upperLimit, directionality); - -} +// void SearchVertex_CallThread(int iStart, int iIte, int nhits, int tolerance, bool likelihood, double lowerLimit, double upperLimit, int directionality) +// { +// LEAF::GetME()->SearchVertex_thread(iStart, iIte, nhits, tolerance, likelihood, lowerLimit, upperLimit, directionality); +// } /*****************************************************************************************************/ -LEAF::LEAF() { - +LEAF::LEAF() +{ myFitter = this; - fThread = N_THREAD; } -LEAF::~LEAF() { - - //fTrueVtxPosDouble.clear(); +LEAF::~LEAF() +{ delete fRand; - } -LEAF* LEAF::GetME() { - - if ( myFitter ) return myFitter; +LEAF *LEAF::GetME() +{ + if (myFitter) + return myFitter; myFitter = new LEAF(); return myFitter; } -void LEAF::DeleteME() { - - if ( myFitter ) delete myFitter; - +void LEAF::DeleteME() +{ + if (myFitter) delete myFitter; } /*****************************************************************************************************/ -void LEAF::Initialize(const Geometry* lGeometry) { - - fGeometry = lGeometry; - - // Get ID PMT info list - fPMTList = fGeometry->GetPMTList(); - - // Set DarkNoise: - fDarkRate_ns[NormalPMT] = fGeometry->pmt_dark_rate[HKAA::kIDPMT_BnL]; - fDarkRate_ns[MiniPMT] = fGeometry->pmt_dark_rate[HKAA::kIDPMT_3inch]; - - // Set Tank size: - fTankRadius = fGeometry->detector_radius; - fTankHeight = fGeometry->detector_length; - fTankHalfHeight = fTankHeight / 2.; - - this->Init(); +void LEAF::Initialize(const Geometry *lGeometry) +{ + InitConfig(lGeometry); + InitInputs(lGeometry); + LoadSplines(); } -/*****************************************************************************************************/ -void LEAF::Init() { +//Creates the candidates list after Coarse grid search, refines the fit, then outputs the best candidate +std::vector LEAF::FitDirection(const std::vector& fixedVertexPosition, FitterOutput &fOutput, int nhits, bool searchPrior) +{ + //* Init + TStopwatch timer; + timer.Reset(); + timer.Start(); + double *tLimits = new double[4]; + int verbose = VERBOSE; + for (int i = 0; i < 4; i++) tLimits[i] = 2 * fSearchVtxStep; + + //* Fit first step : find candidates for the direction + std::vector candidates; + if(searchPrior) + { + //* Only one candidate, chosen wisely by computing a mean hit direction for hits within a specified residual time window + std::vector dirPriorCart = FitDirectionQuick(fixedVertexPosition, fOutput); + std::vector dirPrior = CartesianToPolarNorm(dirPriorCart); + DirectionCandidate priorCandidate; + priorCandidate.theta = dirPrior[0]; + priorCandidate.phi = dirPrior[1]; + candidates.push_back(priorCandidate); + } + else candidates = SearchDirection(fixedVertexPosition, DirTolerance); //* multiple candidates, spread uniformly - // Parameters - fStepByStep = false; // Step by step mode was not test and is not supported by multithreading - fUseDirectionality = false; - fHighEnergy = false; - - // If true, we apply same cuts for timing as B&L PMT when searching vertex. - // Otherwise, we do not apply them and use all the hits. - // The latter is particularly useful when directionality is added, as it can kill DR hits. - fLimit_mPMT = true; - - fSTimePDFLimitsQueueNegative = -3.; //-5;//-3;//-1.5;//-3;//-10; - fSTimePDFLimitsQueuePositive = 4.; //10;//4;//3;//4;//500; - fSTimePDFLimitsQueueNegative_fullTimeWindow = 0.; //-1.5;//-3;//-10; - fSTimePDFLimitsQueuePositive_fullTimeWindow = 0.; //3;//4;//500; - fTimeWindowSizeFull = 1500; - fAveraging = 20;//Number of points we used to average the vertex on. - - fMinimizeLimitsNegative = -700; - fMinimizeLimitsPositive = 1000; - - fSearchVtxStep = 300;//in cm, the step size for coarse grid search - fSearchVtxTolerance = 60;//Number of candidate vertex that are kept after coarse grid search - fHitTimeLimitsNegative = -5;//-10;//-5;//-5 - fHitTimeLimitsPositive = 7;//15;//7;//7 - - fIntegrationTimeWindow = 50;//in ns - - // Compute light speed: - float fCVacuum = 3e8*1e2 / 1e9;//speed of light, in centimeter per ns. - float fNIndex = 1.373;//1.385;//1.373;//refraction index of water - fLightSpeed = fCVacuum/fNIndex; - - /* - fTrueVtxPos.clear(); - fTrueVtxPosDouble.clear(); - fTrueVtxPos.resize(5,0.); - fTrueVtxPosDouble.push_back(fTrueVtxPos); - */ - - fPDFNorm_fullTimeWindow = 0.; - - //fPositionSize = 0; - - // Get random generator - fRand = new TRandom3(); - - this->LoadSplines(); + // std::cout << "nb direction candidates : " << candidates.size() << std::endl; - // Make position lists - this->MakePositionList(); -} + fOutput.Dir_Search_ComputeTime = timer.RealTime(); + TStopwatch timer2; + timer2.Reset(); + timer2.Start(); + //? To be removed for real data where true dir is not known, no impact on the fit + fOutput.stepOneContainsTrueDir = (int)ContainsTrueDir(&candidates, fTrueDir); //* a quick check to see if the step one works (one of the candidate is near the true direction) + if (verbose >= 2) + { + for (unsigned int j = 0; j < candidates.size(); j++) std::cout << "Candidate theta = " << candidates[j].theta * TMath::RadToDeg() << " deg, phi = " << candidates[j].phi * TMath::RadToDeg() << " deg, NLL = " << candidates[j].DNLL << std::endl; + } -void LEAF::LoadSplines() { + DirectionCandidate finalCandidate = MinimizeDirection(fixedVertexPosition, candidates, nhits, tLimits, verbose); - TFile *fSplines, *fSplines2; - if(fHighEnergy){ - fSplines = new TFile("${LEAFDIR}/inputs/timePDF_HE.root","read");//To generate with my code ProduceWSPlots.c - fSplines2 = fSplines; - } - else{ - fSplines = new TFile("${LEAFDIR}/inputs/timePDF_DRnew_Large.root","read");//To generate with my code ProduceWSPlots.c - fSplines2 = new TFile("${LEAFDIR}/inputs/timePDF_Directionality_DRnew.root","read");//To generate with my code ProduceWSPlots.c - } - - // Prevent TGraph2D to be append to TFile (this is needed as we are doing multiple copy of TGraph2D) - // For a strange reason sometimes the TGraph2D is see as a TH1 and stored in an open TFile - TH1::AddDirectory(false); - - for(int pmtType=0;pmtTypeGet(Form("splineExpoQueue%d_%d",0,pmtType)); - fSplineTimePDFDarkRate[pmtType] = (TSpline3*) fSplines->Get(Form("splineDR%d_%d",0,pmtType)); - - fSTimePDFLimitsQueueNegative_fullTimeWindow = fSplineTimePDFQueue[pmtType]->GetXmin(); - fSTimePDFLimitsQueuePositive_fullTimeWindow = fSplineTimePDFQueue[pmtType]->GetXmax(); - //std::cout<<"Min="<Get(Form("fDistResponsePMT_pmtType%d",pmtType)); - } -} + fOutput.Dir = PolarToCartesianNorm({finalCandidate.theta, finalCandidate.phi}); + fOutput.DNLL = finalCandidate.DNLL; -/* -void LEAF::SetTrueVertexInfo(std::vector vtx, double time) { + fOutput.Dir_Minimize_ComputeTime = timer2.RealTime(); - fTrueVtxPosDouble.clear(); - fTrueVtxPos.clear(); - - fTrueVtxPos.resize(5,0.); - - fTrueVtxPos[0] = vtx[0]; - fTrueVtxPos[1] = vtx[1]; - fTrueVtxPos[2] = vtx[2]; - fTrueVtxPos[3] = time; - fTrueVtxPos[4] = 0.; - - fTrueVtxPosDouble.push_back(fTrueVtxPos); -}*/ -/*****************************************************************************************************/ -// Spline Integral functions: -double LEAF::SplineIntegral(TSpline3 * s,double start,double end,double stepSize){ - double integral=0; - for(double i=start;iEval(i+stepSize/2); - } - return integral; -} -double LEAF::SplineIntegralAndSubstract(TSpline3 * s0,TSpline3 * s1,double start,double end,double stepSize){ - double integral=0; - for(double i=start;iEval(i+stepSize/2)-s1->Eval(i+stepSize/2)); - } - return integral; + return fOutput.Dir; } -double LEAF::SplineIntegralExpo(TSpline3 * s,double start,double end,double sigma,double stepSize){ - double integral=0; - for(double i=start;iEval(i+stepSize/2)*TMath::Gaus(i+stepSize/2,0,sigma); - } - return integral; -} +std::vector LEAF::FitDirectionQuick(const std::vector& fixedVertexPosition, FitterOutput &fOutput) +{ + // Quick fit implementation -/*****************************************************************************************************/ - -void LEAF::VectorVertexPMT(std::vector vertex, int iPMT, double* dAngles) { - - // Guillaume 2020/05/20: - // mPMT referencial is computed once for all PMT in LoadPMTInfo() and MakeMPMTReferencial(iPMT) - // Guillaume 2020/11/20: - // Reference is moved to outside LEAF (HKManager or hk-AstroAnalysis) - - PMTInfo lPMTInfo = (*fPMTList)[iPMT]; - int iPMTTop = lPMTInfo.mPMT_RefTube; - PMTInfo lPMTInfoTop = (*fPMTList)[iPMTTop]; - - //5. Now we have our referential, we should just calculate the angles of the PMT to vertex position vector in this referential. - //a. calculate the PMT to vertex position vector. - - double dVtx_PMTRef[3]; - - dVtx_PMTRef[0] = vertex[0] - lPMTInfoTop.Position[0]; - dVtx_PMTRef[1] = vertex[1] - lPMTInfoTop.Position[1]; - dVtx_PMTRef[2] = vertex[2] - lPMTInfoTop.Position[2]; - - double dLengthVtx = Astro_GetLength(dVtx_PMTRef); - GeoTools::Normalize(dVtx_PMTRef); - - if( VERBOSE >= 3 ){ - std::cout << "Vertex position = " << vertex[0] << ", " << vertex[1] << ", " << vertex[2] << std::endl; - std::cout << "Vertex position from PMT = " << dVtx_PMTRef[0] << ", " << dVtx_PMTRef[1] << ", " << dVtx_PMTRef[2] << std::endl; - } - - //b. Then extract Theta and Phi: - double dCosTheta= Astro_GetScalarProd(dVtx_PMTRef,lPMTInfo.mPMT_RefZ); - double dTheta = TMath::ACos(dCosTheta); - - double dPhi = 0.; - - if( (*fPMTList)[iPMT].mPMT_TubeNum == HKAA::kmPMT_TopID ) { - //Phi is not defined in that case.. - } - else { - //We know x=cosPhi x sinTheta and y=sinPhi x sinTheta - double dX = Astro_GetScalarProd(dVtx_PMTRef,lPMTInfo.mPMT_RefX); - double dY = Astro_GetScalarProd(dVtx_PMTRef,lPMTInfo.mPMT_RefY); - double dTanPhi = dY/dX; - dPhi = TMath::ATan(dTanPhi); - - //tan is symetric from -pi/2 to +pi/2 - if( dX == 0 ){ - if( dY < 0 ) dPhi = -TMath::Pi()/2.; - else dPhi = TMath::Pi()/2.; - } - if( dX < 0 ) dPhi += TMath::Pi(); //With this, angle become defines between -pi/2 to -pi/2 +2pi. - - //We wish to bring this from 0 to 2pi: - if( dPhi < 0 ) dPhi += 2*TMath::Pi(); - - //Actually, we have a symmetry in Phi between [0,pi] and [pi,2pi]. So, we will just define Phi in [0,pi] modulo pi - if( dPhi > TMath::Pi() )dPhi = TMath::Pi() - (dPhi - TMath::Pi()); - } - - dAngles[0] = dPhi * 180./TMath::Pi(); - dAngles[1] = dTheta * 180./TMath::Pi(); - dAngles[2] = dLengthVtx; - dAngles[3] = 1; //calculateWeight(dLengthVtx,PMTradius[pmtType],Theta,verbose); - - if( VERBOSE >= 3 ){ - std::cout << "Angles Phi = " << dAngles[0] << ", Theta = " << dAngles[1] << std::endl; - } - -} -double LEAF::FindDirectionTheta(std::vector vertex,int tubeNumber, int verbose) { - - //clock_t timeStart=clock(); - - PMTInfo lPMTInfo = (*fPMTList)[tubeNumber]; - - int pmt_number_in_mpmt = lPMTInfo.mPMT_TubeNum; - double particleRelativePMTpos[3]; - for(int j=0;j<3;j++) particleRelativePMTpos[j] = lPMTInfo.Position[j] - vertex[j]; - if(verbose == 3){ - std::cout<<"Original PMT number = "<MakeEventInfo(lowerLimit,upperLimit); - - //First, substract the DR from the PDF to keep only the signal: - double DR=fSplineTimePDFDarkRate[pmtType]->Eval(residual); - proba-=DR; - //Then, scale signal so that signal integral / DR integral = signalOverNoise - //To do so, we should scale signal so that integral = DR integral * signalOverNoise - //And of course, to rescale, we should divide by the signal integral - double Factor = fEventInfo[pmtType].NoiseIntegral * fEventInfo[pmtType].SignaloverNoise / fEventInfo[pmtType].SignalIntegral; - proba*=Factor; - //And add again the DR - proba+=DR; -#ifdef VERBOSE_NLL - if(VERBOSE>=3) std::cout<<"proba after scaling="<= upperLimit ){ - proba=fSplineTimePDFQueue[pmtType]->Eval(upperLimit); -#ifdef VERBOSE_NLL - if(VERBOSE>=3 && pmtType==1) std::cout<<"PMT type = "<< pmtType <<", Upper limit = "<-1e-1) proba=0;//Since spline sometimes slightly goes below 0 due to interpolation - else { - std::cout <<"Error in "<< residual << "ns where proba = " << proba << std::endl; - return 0; - } - } - - if(proba==0) proba=1e-20; - NLL += -TMath::Log(proba); - } - - //std::cout << " FindNLL_Likelihood loop: " << timer.RealTime() << " " << nhits << std::endl; - if(directionality!=0){ - //double NLLdirBaye=0; - double NLLdir=0; - - NLLdir=this->FindNLLDirectionality(vertexPosition, nhits, VERBOSE,lowerLimit,upperLimit); - -#ifdef VERBOSE_NLL - if(VERBOSE>=2) std::cout<<"NLL = "< vertexPosition, int nhits, double /*lowerLimit*/, double /*upperLimit*/, bool /*killEdges*/, bool /*scaleDR*/, int /*directionality*/){ - double NLL = 0; - - //std::cout << " Find NLL " << nhits << std::endl; - //std::cout << " First Hit " << fHitCollection->At(0).PMT << " " << vertexPosition.size() << std::endl; - for(int ihit = 0; ihit < nhits; ihit++){ - //Hit lHit = fHitInfo[ihit]; - Hit lHit = fHitCollection->At(ihit); - - //std::cout << " NLL Hit " << ihit << " " << lHit.PMT << std::endl; - int iPMT = lHit.PMT; - - double hitTime = (fTimeCorrection + lHit.T) / TimeDelta::ns; - PMTInfo lPMTInfo = (*fPMTList)[iPMT]; - - int pmtType = Astro_GetPMTType(iPMT); - double distance = Astro_GetDistance(lPMTInfo.Position,vertexPosition); - - double tof = distance / fLightSpeed; - double residual = hitTime - tof - vertexPosition[3]; -#ifdef KILLHALF - if(pmtType==1){ - double t=fRand->Uniform(0,1); - if(t>0.5) continue; - } -#endif + // for(int j = 0; j < 3 ; j++) + // { + // if(pmtPosition[j] > maxPos[j]) maxPos[j] = pmtPosition[j]; + // if(pmtPosition[j] < minPos[j]) minPos[j] = pmtPosition[j]; + // } + } - bool bCondition = (residual > fHitTimeLimitsNegative && residual < fHitTimeLimitsPositive) || (pmtType == 1 && !fLimit_mPMT); - - //std::cout << " HIT " << ihit << " Has " << bCondition << " " << pmtType << " " << fLimit_mPMT << " " << residual << " ( " << hitTime << " - " << tof << " - " << vertexPosition[3] << " ) " << distance << " / " << fLightSpeed << " " << lPMTInfo.Position[0] << " "<< lPMTInfo.Position[1] << " "<< lPMTInfo.Position[2] << " " << iPMT << std::endl; - - if( bCondition ){ - NLL++; -#ifdef VERBOSE_NLL - if(VERBOSE>=3){ - std::cout<<"Residual="< meanPos(3, 0.0); + for(long unsigned int i = 0; i < pmtPositions.size(); i++) + { + if(residuals[i] < 10 && residuals[i] > -10) + { + for(int j = 0; j < 3; j++) meanPos[j] += pmtPositions[i][j] * pmtCharges[i]; + sumCharge += pmtCharges[i]; } } - - //std::cout << " NLL Final" << NLL << std::endl; - - //std::cout<<"NLL="<MakeEventInfo(lowerLimit,upperLimit); - - //First, substract the DR from the PDF to keep only the signal: - double DR=fSplineTimePDFDarkRate[pmtType]->Eval(residual); - proba-=DR; - //Then, scale signal so that signal integral / DR integral = signalOverNoise - //To do so, we should scale signal so that integral = DR integral * signalOverNoise - //And of course, to rescale, we should divide by the signal integral - double Factor = fEventInfo[pmtType].NoiseIntegral * fEventInfo[pmtType].SignaloverNoise / fEventInfo[pmtType].SignalIntegral; - proba*=Factor; - //And add again the DR - proba+=DR; - if(VERBOSE>=3) std::cout<<"proba after scaling="<= upperLimit){ - proba=fSplineTimePDFQueue[pmtType]->Eval(upperLimit); -#ifdef VERBOSE_NLL - if(VERBOSE>=3 && pmtType==1) std::cout<<"PMT type = "<< pmtType <<", Upper limit = "<-1e-1) proba=0;//Since spline sometimes slightly goes below 0 due to interpolation - else { - std::cout <<"Error in "<< residual << "ns where proba = " << proba << std::endl; - return 0; - } - } - if(proba==0) proba=1e-20; - NLL += -TMath::Log(proba); - } - } - std::cout<<"FindNLL: Time it took for the loop = "<< timer.RealTime() << std::endl; - //cout<<"NLL="<=2) cout<<"NLL="<FindNLLDirectionality(vertexPosition, nhits, verbose,lowerLimit,upperLimit); - } - else{ - //NLLdir=findNLLDirectionalityBayes(vertexPosition, nhits, verbose,fHitTimeLimitsNegative,fHitTimeLimitsPositive); - NLLdir2=this->FindNLLDirectionality(vertexPosition, nhits, verbose,fHitTimeLimitsNegative,fHitTimeLimitsPositive); - } - if(VERBOSE>=2) std::cout<<"NLL = "<SetErrorDef(1); + + if (verbose < 2) minuit->mnexcm("SET NOWarnings", 0, 0, err); + + arglist[0] = 2; + minuit->mnexcm("SET STR", arglist, 1, err); + minimizer->SetFCN(MinuitDirNLL); + + for (int icand = 0; icand < (int)candidates.size(); icand++) + { + // Set initial parameters and bounds + minimizer->SetParameter(0, "theta", candidates[icand].theta, 2*M_PI/180, 0.0, TMath::Pi()); + minimizer->SetParameter(1, "phi", candidates[icand].phi, 2*M_PI/180, -TMath::Pi(), TMath::Pi()); + minimizer->SetParameter(2, "nhits", nhits, nhits, nhits - 1, nhits + 1); + minimizer->SetParameter(3, "vertex0", fixedVertexPosition[0], 5e-1, -limits[0], limits[0]); + minimizer->SetParameter(4, "vertex1", fixedVertexPosition[1], 5e-1, -limits[1], limits[1]); + minimizer->SetParameter(5, "vertex2", fixedVertexPosition[2], 5e-1, -limits[2], limits[2]); + minimizer->SetParameter(6, "vertex3", fixedVertexPosition[3], 5e-1 / fLightSpeed, -limits[3] / fLightSpeed, limits[3] / fLightSpeed); + + minimizer->FixParameter(2); + minimizer->FixParameter(3); + minimizer->FixParameter(4); + minimizer->FixParameter(5); + minimizer->FixParameter(6); - return NLL; -} + // Store fixed vertex position and number of hits in member variables + // fFixedVertexPosition = fixedVertexPosition; + // fNHits = nhits; -double LEAF::FindNLLDirectionality(std::vector vVtxPos, int nhits, int verbose, double /*lowerLimit*/, double /*upperLimit*/) { + // Minimize + arglist[0] = 1e9; // Maximum number of function evaluations + arglist[1] = 1e-1; // Tolerance + minimizer->ExecuteCommand("MIGRAD", arglist, 2 ); - double NLL = 0; - std::vector vDirection(2,0.);//Return phi and theta. - - // Interpolate isn't compatible with multi-thread, need to use mutex which lead to long deadtime. - // Copy TGraph2D - mtx.lock(); - static thread_local TGraph2D gPMTDirectionality_2D_local_0 = TGraph2D(*gPMTDirectionality_2D[MiniPMT][0]); - static thread_local TGraph2D gPMTDirectionality_2D_local_1 = TGraph2D(*gPMTDirectionality_2D[MiniPMT][1]); - static thread_local TGraph2D gPMTDirectionality_2D_local_2 = TGraph2D(*gPMTDirectionality_2D[MiniPMT][2]); - mtx.unlock(); - - TGraph2D* tDirectionality[3]; - tDirectionality[0] = &gPMTDirectionality_2D_local_0; - tDirectionality[1] = &gPMTDirectionality_2D_local_1; - tDirectionality[2] = &gPMTDirectionality_2D_local_2; - - - for(int ihit = 0; ihit < nhits; ihit++){ - //Hit lHit = fHitInfo[ihit]; - Hit lHit = fHitCollection->At(ihit); - - int iPMT = lHit.PMT; - - PMTInfo lPMTInfo = (*fPMTList)[iPMT]; - int pmtType = Astro_GetPMTType(iPMT); + // Retrieve results + std::vector optimizedParameters(2); + candidates[icand].theta = minimizer->GetParameter(0); // theta + candidates[icand].phi = minimizer->GetParameter(1); // phi - if(pmtType==0) continue; - bool condition = true; - - if(condition){ - double vPMTVtx[4]; - this->VectorVertexPMT(vVtxPos,iPMT,vPMTVtx); - - double dPhi = vPMTVtx[0]; - double dTheta = vPMTVtx[1]; - double dDist = vPMTVtx[2]; - - int pmtGroup = lPMTInfo.mPMT_Group; - - double proba = 0; - proba = tDirectionality[pmtGroup]->Interpolate(dPhi,dTheta); - - double dDistCorr = fDistResponsePMT[pmtType]->Eval(dDist); - - proba *= dDistCorr; - - if(proba==0){ - proba=1e-20; - if(verbose) std::cout<<"We are at proba = 0, theta = "<< dTheta <<", proba used = "<< proba <=3){ - //std::cout<<"PMT type="< > LEAF::SearchVertex(int nhits,int tolerance,bool likelihood,double lowerLimit, double upperLimit,int directionality){ - //2. How to set the search? - //double stepSize = 0.5;//in m - +// Coarse search of the vertex. It will search in a cylinder having the radius = tankRadius and a height tankHeight. +// The step between each grid points is given by stepSize (in centimeters). +// We give the list of hit information through the 2D array fHitInfo to make the code faster, instead of using the whole code information +// The code provide as an output not only one vertex, but a list of vertices whose likelihood is higher (Negative Log Likelihood is lower). +// The number of output vertices can be chosen using "tolerance". +std::vector> LEAF::SearchVertex(int nhits, int tolerance, bool likelihood, double lowerLimit, double upperLimit, int directionality) +{ + // 2. How to set the search? + // double stepSize = 0.5;//in m + std::vector tBestReconstructedVertexPosition; - - //double ** reconstructedVertexPosition = new double[4]; - - std::vector tVtxContainer; - - clock_t timeStart=clock(); - - if ( fPositionList.size() == 0 ) std::cout << " Position list not filled " << std::endl; - if(VERBOSE>=3) std::cout<<"Number of hits for coarse grid search = "<FindNLL_Likelihood (hPos.Vtx,nhits,lowerLimit,upperLimit,false,false,directionality); - else hPos.NLL = this->FindNLL_NoLikelihood(hPos.Vtx,nhits,lowerLimit,upperLimit,false,false,directionality); - - hPos.NLL += fRand->Uniform(0,1e-5);//Here only not to have exactly the same value for 2NLL, and be removed from map - - tVtxContainer.push_back(hPos); - - if(VERBOSE>=3) std::cout << "Test Pos: " << hPos.Vtx[3] << " " << hPos.Vtx[0] << " " << hPos.Vtx[1] << " " << hPos.Vtx[2] << " NLL = " << hPos.NLL << " " << likelihood << std::endl; - + + // double ** reconstructedVertexPosition = new double[4]; + + std::vector tVtxContainer; + + clock_t timeStart = clock(); + + if (fPositionList.size() == 0) + std::cout << " Position list not filled " << std::endl; + if (VERBOSE >= 3) + std::cout << "Number of hits for coarse grid search = " << nhits << std::endl; + for (unsigned int iPos = 0; iPos < fPositionList.size(); iPos++) + { + struct FitPosition hPos; + hPos.Vtx = fPositionList[iPos]; + hPos.NLL = 0; + + if (likelihood) + hPos.NLL = Likelihoods::Vertex_Time_NLL(fHitCollection, hPos.Vtx, nhits, lowerLimit, upperLimit, false, false, directionality); + else + hPos.NLL = Likelihoods::Vertex_Score(fHitCollection, hPos.Vtx, nhits, lowerLimit, upperLimit, false, false, directionality); + + hPos.NLL += fRand->Uniform(0, 1e-5); // Here only not to have exactly the same value for 2NLL, and be removed from map + + tVtxContainer.push_back(hPos); + + if (VERBOSE >= 3) + std::cout << "Test Pos: " << hPos.Vtx[3] << " " << hPos.Vtx[0] << " " << hPos.Vtx[1] << " " << hPos.Vtx[2] << " NLL = " << hPos.NLL << " " << likelihood << std::endl; + #ifdef VERBOSE_VTX - //if(verbose && distanceToTrue < 2) std::cout << "Distance to true vertex = " << distanceToTrue << ", Probability = " << NLL < > tReconstructedVertexPosition; - for ( int iPos = 0; iPos < tolerance; iPos++ ) { - std::vector vPos(5,0.); + + clock_t timeEndLoop = clock(); + + std::vector> tReconstructedVertexPosition; + for (int iPos = 0; iPos < tolerance; iPos++) + { + std::vector vPos(5, 0.); struct FitPosition hPos = tVtxContainer[iPos]; - + vPos[0] = hPos.Vtx[0]; vPos[1] = hPos.Vtx[1]; vPos[2] = hPos.Vtx[2]; vPos[3] = hPos.Vtx[3]; vPos[4] = hPos.NLL; - + tReconstructedVertexPosition.push_back(vPos); - + #ifdef VERBOSE_VTX - std::cout<<"NLL = " << hPos.NLL << ", vertex = " << vPos[0]<<", "< > LEAF::SearchVertex_Main( - int nhits,int tolerance,bool likelihood,double lowerLimit, double upperLimit,int directionality){ - - +std::vector LEAF::SearchVertex_Main(int nhits, int tolerance, bool likelihood, double lowerLimit, double upperLimit, int directionality) +{ +// 1) Launch threads just like in SearchVertex_Main. std::vector lThreadList; - std::vector< std::vector > lOutputFinal; - - int iCand_Step = (int) fPositionList.size() / fThread; - - if ( iCand_Step < 1 ) iCand_Step = 1; - - for ( int iStart = 0; iStart < (int) fPositionList.size(); iStart += iCand_Step ) { - - std::thread tThrd( SearchVertex_CallThread, iStart,iCand_Step, - nhits, tolerance, likelihood, lowerLimit, upperLimit, directionality); - + std::vector> lOutputFinal; + TStopwatch timer; + timer.Reset(); + timer.Start(); + + int iCand_Step = (int)fPositionList.size()/fThread; + if(iCand_Step < 1) iCand_Step = 1; + + for(int iStart=0; iStart<(int)fPositionList.size(); iStart+=iCand_Step) + { + std::thread tThrd( + &LEAF::SearchVertex_thread, + LEAF::GetME(), + iStart, iCand_Step, + nhits, tolerance, + likelihood, lowerLimit, upperLimit, directionality + ); lThreadList.push_back(std::move(tThrd)); - } - - - for ( unsigned int iIdx=0; iIdx < lThreadList.size(); iIdx++ ) { - // Wait thread - - if ( lThreadList[iIdx].joinable() ) - lThreadList[iIdx].join(); - } - - // Check if mutex is already lock, if not -> lock - // We don't want simultaneous modification of the common output + + // std::cout << "time after vtx search threads in SearchVertex_Main : " << timer.RealTime() << std::endl; + + for(auto &th: lThreadList) if(th.joinable()) th.join(); + + // 2) collect partial results from fThreadOutput. mtx.lock(); - lOutputFinal = fThreadOutput; + lOutputFinal = fThreadOutput; // each row: [x, y, z, t, NLL] fThreadOutput.clear(); mtx.unlock(); - - std::sort(lOutputFinal.begin(), lOutputFinal.end(), SortOutputVector); - while((int) lOutputFinal.size() > tolerance) { - lOutputFinal.pop_back(); - } - - //std::cout << " SearchVtx main " << tolerance << " " << lOutputFinal.size() << std::endl; - return lOutputFinal; -} + // 3) Sort by NLL and keep 'tolerance' number of candidates + std::sort(lOutputFinal.begin(), lOutputFinal.end(), SortOutputVector); + while((int)lOutputFinal.size() > tolerance) lOutputFinal.pop_back(); + + //4) Now build a vector of VtxCandidate. + std::vector finalList; + finalList.reserve(lOutputFinal.size()); + + for(auto &cand : lOutputFinal) + { + // cand = [ x, y, z, t, NLL ] + // double vx = cand[0]; + // double vy = cand[1]; + // double vz = cand[2]; + // double vt = cand[3]; + // double nll= cand[4]; + + // Build a 4-vector. + // std::vector vtx(4); + // vtx[0] = vx; + // vtx[1] = vy; + // vtx[2] = vz; + // vtx[3] = vt; + + // 5) Per-candidate SNR. + VtxCandidate cOut = CreateCandidate(cand, lowerLimit, upperLimit, false); + cOut.NLL = cand[4]; // store the NLL from the coarse search. + finalList.push_back(cOut); + } + + // std::cout << "time at end of vertex search in SearchVertex_Main : " << timer.RealTime() << std::endl; + return finalList; +} + +// Coarse search of the vertex. It will search in a cylinder having the radius = tankRadius and a height tankHeight. +// The step between each grid points is given by stepSize (in centimeters). +// We give the list of hit information through the 2D array fHitInfo to make the code faster, instead of using the whole code information +// The code provide as an output not only one vertex, but a list of vertices whose likelihood is higher (Negative Log Likelihood is lower). +// The number of output vertices can be chosen using "tolerance". +void LEAF::SearchVertex_thread(int iStart, int iIte, int nhits, int tolerance, bool likelihood, double lowerLimit, double upperLimit, int directionality) +{ + // 2. How to set the search? + // double stepSize = 0.5;//in m + + std::vector tVtxContainer; + + TStopwatch timer; + timer.Reset(); + timer.Start(); + + if (fPositionList.size() == 0) + std::cout << " Position list not filled " << std::endl; + if (VERBOSE >= 3) + std::cout << "Number of hits for coarse grid search = " << nhits << std::endl; + int iEnd = (iIte + iStart); + if (iEnd > (int)fPositionList.size()) + iEnd = fPositionList.size(); -//Coarse search of the vertex. It will search in a cylinder having the radius = tankRadius and a height tankHeight. -//The step between each grid points is given by stepSize (in centimeters). -//We give the list of hit information through the 2D array fHitInfo to make the code faster, instead of using the whole code information -//The code provide as an output not only one vertex, but a list of vertices whose likelihood is higher (Negative Log Likelihood is lower). -//The number of output vertices can be chosen using "tolerance". -void LEAF::SearchVertex_thread( - int iStart, int iIte, - int nhits,int tolerance,bool likelihood,double lowerLimit, double upperLimit,int directionality){ - //2. How to set the search? - //double stepSize = 0.5;//in m - - std::vector tVtxContainer; - - if ( fPositionList.size() == 0 ) std::cout << " Position list not filled " << std::endl; - if(VERBOSE>=3) std::cout<<"Number of hits for coarse grid search = "< (int) fPositionList.size() ) iEnd = fPositionList.size(); - - for ( int iPos = iStart; iPos < iEnd; iPos++ ) { - - struct FitPosition hPos; - hPos.Vtx = fPositionList[iPos]; - hPos.NLL = 0; - //std::cout << " Thread2 " << iPos << " " << iPos-iStart << " " << iEnd << std::endl; - - if ( likelihood ) hPos.NLL = this->FindNLL_Likelihood (hPos.Vtx,nhits,lowerLimit,upperLimit,false,false,directionality); - else hPos.NLL = this->FindNLL_NoLikelihood(hPos.Vtx,nhits,lowerLimit,upperLimit,false,false,directionality); - - //std::cout << " Thread3 " << iPos << " " << iPos-iStart << " " << iEnd << std::endl; - //hPos.NLL += fRand->Uniform(0,1e-5);//Here only not to have exactly the same value for 2NLL, and be removed from map - - tVtxContainer.push_back(hPos); - } - - std::sort(tVtxContainer.begin(),tVtxContainer.end(),SortingNLL()); - while((int) tVtxContainer.size() > tolerance) { - tVtxContainer.pop_back(); + for (int iPos = iStart; iPos < iEnd; iPos++) + { + + struct FitPosition hPos; + hPos.Vtx = fPositionList[iPos]; + hPos.NLL = 0; + // std::cout << " Thread2 " << iPos << " " << iPos-iStart << " " << iEnd << std::endl; + + if (likelihood) hPos.NLL = Likelihoods::Vertex_Time_NLL(fHitCollection, hPos.Vtx, nhits, lowerLimit, upperLimit, false, false, directionality); + else hPos.NLL = Likelihoods::Vertex_Score(fHitCollection, hPos.Vtx, nhits, lowerLimit, upperLimit, false, false, directionality); + + tVtxContainer.push_back(hPos); } - + + std::sort(tVtxContainer.begin(), tVtxContainer.end(), Likelihoods::SortingNLL()); + while ((int)tVtxContainer.size() > tolerance) tVtxContainer.pop_back(); + // Check if mutex is already lock, if not -> lock // We don't want simultaneous modification of the common output mtx.lock(); - //std::vector< std::vector > tReconstructedVertexPosition; - for ( unsigned int iPos = 0; iPos < tVtxContainer.size(); iPos++ ) { - std::vector vPos(5,0.); + // std::vector< std::vector > tReconstructedVertexPosition; + for (unsigned int iPos = 0; iPos < tVtxContainer.size(); iPos++) + { + std::vector vPos(5, 0.); struct FitPosition hPos = tVtxContainer[iPos]; - + vPos[0] = hPos.Vtx[0]; vPos[1] = hPos.Vtx[1]; vPos[2] = hPos.Vtx[2]; @@ -987,80 +523,96 @@ void LEAF::SearchVertex_thread( fThreadOutput.push_back(vPos); } mtx.unlock(); - -} + // std::cout << "thread took : " << timer.RealTime() << std::endl; +} -//Fine grid search of the vertex. It will search in a box this time, and not a cylinder as searchVertex is doing. It is useful when we have one or several first candidate vertices. Therefore, this code is generally run after serachVertex. -//The box is centered around the initial vertex provided. -//The box size is 2 x limits in each direction, centered on the initial vertex provided. -//nCandidate provide the size of the list of inital vertices, while tolerance provide the size of the output list. -//Other informations are the same as searchVertex. -std::vector< std::vector > LEAF::SearchVertexFine(std::vector< std::vector > initialVertex,double * limits, double stepSize,int nhits,int nCandidates,int tolerance,int verbose,bool likelihood,bool average,double lowerLimit, double upperLimit, int directionality){ - - if(verbose) std::cout<<"Start fine search, use directionality? "<> LEAF::SearchVertexFine(std::vector> initialVertex, double *limits, double stepSize, int nhits, int nCandidates, int tolerance, int verbose, bool likelihood, bool average, double lowerLimit, double upperLimit, int directionality) +{ + if (verbose) std::cout << "Start fine search, use directionality? " << directionality << std::endl; + // 2. How to set the search? + // double stepSize = 0.5;//in m + // double * reconstructedVertexPosition = new double[4]; std::vector tBestReconstructedVertexPosition; - double minNLL=1e15; - clock_t timeStart=clock(); - - std::map< double, std::vector > vertexContainer; - - for(int icand=0;icand=2) std::cout << "Time = " << time << "ns" << std::endl; - for(double x = initialVertex[icand][VTX_X] - limits[1]; x <= initialVertex[icand][VTX_X] + limits[1];x+=stepSize){ - if(VERBOSE>=2) std::cout << "x = " << x << "m" << std::endl; - for(double y = initialVertex[icand][VTX_Y] - limits[2]; y <= initialVertex[icand][VTX_Y] + limits[2];y+=stepSize){ - if(VERBOSE>=2) std::cout << "y = " << y << "m" << std::endl; - for(double z = initialVertex[icand][VTX_Z] - limits[3]; z <= initialVertex[icand][VTX_Z] + limits[3];z+=stepSize){ - if(VERBOSE>=2) std::cout << "z = " << z << "m" << std::endl; - - std::vector tVtxPosition(4,0.);//In centimeters - tVtxPosition[0] = x;//radius*TMath::Cos(angle); - tVtxPosition[1] = y;//radius*TMath::Sin(angle); - tVtxPosition[2] = z;//height; + double minNLL = 1e15; + clock_t timeStart = clock(); + + std::map> vertexContainer; + + for (int icand = 0; icand < nCandidates; icand++) + { + if (verbose) + std::cout << "Candidate #" << icand << std::endl; + for (double time = initialVertex[icand][VTX_T] - limits[0] / fLightSpeed; time <= initialVertex[icand][VTX_T] + limits[0] / fLightSpeed; time += (stepSize / fLightSpeed)) + { + if (VERBOSE >= 2) + std::cout << "Time = " << time << "ns" << std::endl; + for (double x = initialVertex[icand][VTX_X] - limits[1]; x <= initialVertex[icand][VTX_X] + limits[1]; x += stepSize) + { + if (VERBOSE >= 2) + std::cout << "x = " << x << "m" << std::endl; + for (double y = initialVertex[icand][VTX_Y] - limits[2]; y <= initialVertex[icand][VTX_Y] + limits[2]; y += stepSize) + { + if (VERBOSE >= 2) + std::cout << "y = " << y << "m" << std::endl; + for (double z = initialVertex[icand][VTX_Z] - limits[3]; z <= initialVertex[icand][VTX_Z] + limits[3]; z += stepSize) + { + if (VERBOSE >= 2) + std::cout << "z = " << z << "m" << std::endl; + + std::vector tVtxPosition(4, 0.); // In centimeters + tVtxPosition[0] = x; // radius*TMath::Cos(angle); + tVtxPosition[1] = y; // radius*TMath::Sin(angle); + tVtxPosition[2] = z; // height; tVtxPosition[3] = time; - double NLL=0; - if(average){ - //double vNLL[fAveraging]; - for(int irand=0;irandUniform(0,stepSize); - double phi = fRand->Uniform(0,2*TMath::Pi()); - double theta = fRand->Uniform(0,TMath::Pi()); - double timeRand = 1.5*fRand->Uniform(-stepSize/fLightSpeed,stepSize/fLightSpeed); - std::vector randVertexPosition(4,0.);//In centimeters - randVertexPosition[0] = tVtxPosition[0] + radius*TMath::Sin(theta)*TMath::Cos(phi); - randVertexPosition[1] = tVtxPosition[1] + radius*TMath::Sin(theta)*TMath::Sin(phi); - randVertexPosition[2] = tVtxPosition[2] + radius*TMath::Cos(theta); - randVertexPosition[3] = tVtxPosition[3] + timeRand;//rand->Uniform(-stepSize/fLightSpeed/2,stepSize/fLightSpeed/2); - //Evaluate its NLL - double nll=this->FindNLL(randVertexPosition,nhits,likelihood,verbose,lowerLimit,upperLimit,false,false,directionality); - NLL+=nll; - //vNLL[irand]=nll; + double NLL = 0; + if (average) + { + // double vNLL[fAveraging]; + for (int irand = 0; irand < fAveraging; irand++) + { + // Determine the position around the vertex + double radius = 1.5 * fRand->Uniform(0, stepSize); + double phi = fRand->Uniform(0, 2 * TMath::Pi()); + double theta = fRand->Uniform(0, TMath::Pi()); + double timeRand = 1.5 * fRand->Uniform(-stepSize / fLightSpeed, stepSize / fLightSpeed); + std::vector randVertexPosition(4, 0.); // In centimeters + randVertexPosition[0] = tVtxPosition[0] + radius * TMath::Sin(theta) * TMath::Cos(phi); + randVertexPosition[1] = tVtxPosition[1] + radius * TMath::Sin(theta) * TMath::Sin(phi); + randVertexPosition[2] = tVtxPosition[2] + radius * TMath::Cos(theta); + randVertexPosition[3] = tVtxPosition[3] + timeRand; // rand->Uniform(-stepSize/fLightSpeed/2,stepSize/fLightSpeed/2); + // Evaluate its NLL + double nll = Likelihoods::FindNLL(fHitCollection, randVertexPosition, nhits, likelihood, verbose, lowerLimit, upperLimit, false, false, directionality); + NLL += nll; + // vNLL[irand]=nll; } - - if(fAveraging!=0) NLL/=fAveraging; + + if (fAveraging != 0) + NLL /= fAveraging; } - else NLL=this->FindNLL(tVtxPosition,nhits,likelihood,verbose,lowerLimit,upperLimit,false,false,directionality); + else + NLL = Likelihoods::FindNLL(fHitCollection, tVtxPosition, nhits, likelihood, verbose, lowerLimit, upperLimit, false, false, directionality); + + std::vector candidateReconstructedVertexPosition(4, 0.); + for (int i = 0; i < 4; i++) + candidateReconstructedVertexPosition[i] = tVtxPosition[i]; - - std::vector candidateReconstructedVertexPosition(4,0.); - for(int i=0;i<4;i++) candidateReconstructedVertexPosition[i] = tVtxPosition[i]; - vertexContainer[NLL] = candidateReconstructedVertexPosition; - if( (int) vertexContainer.size() > tolerance){ + if ((int)vertexContainer.size() > tolerance) + { vertexContainer.erase(--vertexContainer.rbegin().base()); } - - if(NLL> LEAF::MinimizeVertex(std::vector> initialVertex, double *limits, double stepSize, int nhits, int nCandidates, int tolerance, int verbose, bool /*likelihood*/, bool /*average*/, double lowerLimit, double upperLimit, int directionality) +{ + if (VERBOSE >= 2) std::cout << "Minimizer" << std::endl; + // 2. How to set the search? + // double stepSize = 0.5;//in m + // double * reconstructedVertexPosition = new double[4]; std::vector tBestReconstructedVertexPosition; - clock_t timeStart=clock(); + clock_t timeStart = clock(); - std::vector tVtxContainer; + std::vector tVtxContainer; + + if (VERBOSE >= 2) + std::cout << "Minimizer" << std::endl; + TFitter *minimizer = new TFitter(8); // 4=nb de params? + TMinuit *minuit = minimizer->GetMinuit(); - if(VERBOSE>=2) std::cout<<"Minimizer"<GetMinuit(); - double arglist[20]; - int err=0; - //arglist[0]=0; - double p1=verbose-1; - minimizer->ExecuteCommand("SET PRINTOUT",&p1,1);//quiet mode + int err = 0; + // arglist[0]=0; + double p1 = verbose - 1; + minimizer->ExecuteCommand("SET PRINTOUT", &p1, 1); // quiet mode minuit->SetErrorDef(1); - //minuit->SetPrintLevel(-1); // quiet mode - if(VERBOSE<2) { - minuit->mnexcm("SET NOWarnings",0,0,err); - } - arglist[0]=2; - minuit->mnexcm("SET STR",arglist,1,err);//set strategy sets the number of derivative estimated. in the present case (2), this is the highest number of it with a counter-balance: very slow! - minimizer->SetFCN(MinuitLikelihood);//here is function to minimize - - double Mig[2]={1e6,1e0};//maxcalls and tolerance - //double Imp2[1]={10000}; - - minimizer->SetParameter(0,"vertex0",0,stepSize ,-limits[0] ,limits[0] ); - minimizer->SetParameter(1,"vertex1",0,stepSize ,-limits[1] ,limits[1] ); - minimizer->SetParameter(2,"vertex2",0,stepSize ,-limits[2] ,limits[2] ); - minimizer->SetParameter(3,"vertex3",0,stepSize/fLightSpeed,-limits[3]/fLightSpeed ,limits[3]/fLightSpeed ); - - minimizer->SetParameter(4,"PMTConfiguration",0,0,0,0); // Useless now - minimizer->SetParameter(5,"nhits",nhits,nhits,nhits-1,nhits+1); - minimizer->SetParameter(6,"lowerLimit",lowerLimit,5e-1,-7,-2); - minimizer->SetParameter(7,"upperLimit",upperLimit,5e-1,2,10); - //minimizer->SetParameter(8,"scaling factor",1,1,0,1e3); - //minimizer->SetParameter(8,"signal pe",1,1,0,1e3); - minimizer->SetParameter(8,"expoSigma",100,5,0,1e3); - minimizer->SetParameter(9,"directionality",directionality,1,0,2); - minimizer->SetParameter(10,"thread",0,0,0,0); + // minuit->SetPrintLevel(-1); // quiet mode + if (VERBOSE < 2) minuit->mnexcm("SET NOWarnings", 0, 0, err); + arglist[0] = 2; + minuit->mnexcm("SET STR", arglist, 1, err); // set strategy sets the number of derivative estimated. in the present case (2), this is the highest number of it with a counter-balance: very slow! + minimizer->SetFCN(MinuitLikelihood); // here is function to minimize + + double Mig[2] = {1e6, 1e0}; // maxcalls and tolerance + // double Imp2[1]={10000}; + + minimizer->SetParameter(0, "vertex0", 0, stepSize, -limits[0], limits[0]); + minimizer->SetParameter(1, "vertex1", 0, stepSize, -limits[1], limits[1]); + minimizer->SetParameter(2, "vertex2", 0, stepSize, -limits[2], limits[2]); + minimizer->SetParameter(3, "vertex3", 0, stepSize / fLightSpeed, -limits[3] / fLightSpeed, limits[3] / fLightSpeed); + + minimizer->SetParameter(4, "PMTConfiguration", 0, 0, 0, 0); // Useless now + minimizer->SetParameter(5, "nhits", nhits, nhits, nhits - 1, nhits + 1); + minimizer->SetParameter(6, "lowerLimit", lowerLimit, 5e-1, -7, -2); + minimizer->SetParameter(7, "upperLimit", upperLimit, 5e-1, 2, 10); + // minimizer->SetParameter(8,"scaling factor",1,1,0,1e3); + // minimizer->SetParameter(8,"signal pe",1,1,0,1e3); + minimizer->SetParameter(8, "expoSigma", 100, 5, 0, 1e3); + minimizer->SetParameter(9, "directionality", directionality, 1, 0, 2); + minimizer->SetParameter(10, "thread", 0, 0, 0, 0); minimizer->FixParameter(4); minimizer->FixParameter(5); minimizer->FixParameter(6); @@ -1149,31 +702,33 @@ std::vector< std::vector > LEAF::MinimizeVertex(std::vector< std::vector minimizer->FixParameter(9); minimizer->FixParameter(10); - for(int icand=0;icand=1){ - std::cout<<"Candidate #"<= 1) + { + std::cout << "Candidate #" << icand << std::endl; + std::cout << "Initial (x,y,z,t): " << initialVertex[icand][0] + << ", " << initialVertex[icand][1] + << ", " << initialVertex[icand][2] + << ", " << initialVertex[icand][3] << std::endl; } #endif - - minimizer->SetParameter(0,"vertex0",initialVertex[icand][0],stepSize ,initialVertex[icand][0]-limits[0] ,initialVertex[icand][0]+limits[0] ); - minimizer->SetParameter(1,"vertex1",initialVertex[icand][1],stepSize ,initialVertex[icand][1]-limits[1] ,initialVertex[icand][1]+limits[1] ); - minimizer->SetParameter(2,"vertex2",initialVertex[icand][2],stepSize ,initialVertex[icand][2]-limits[2] ,initialVertex[icand][2]+limits[2] ); - minimizer->SetParameter(3,"vertex3",initialVertex[icand][3],stepSize/fLightSpeed ,initialVertex[icand][3]-limits[3]/fLightSpeed ,initialVertex[icand][3]+limits[3]/fLightSpeed ); - - minimizer->SetParameter(4,"PMTConfiguration",0,0,0,0); // Useless now - minimizer->SetParameter(5,"nhits",nhits,nhits,nhits-1,nhits+1); - minimizer->SetParameter(6,"lowerLimit",lowerLimit,5e-1,-7,-2); - minimizer->SetParameter(7,"upperLimit",upperLimit,5e-1,2,10); - minimizer->SetParameter(8,"expoSigma",100,5,0,1e3); - //minimizer->SetParameter(8,"scaling factor",1,1,0,1e3); - //minimizer->SetParameter(8,"signal pe",1,1,0,1e3); - minimizer->SetParameter(9,"directionality",directionality,1,0,2); - minimizer->SetParameter(10,"thread",0,0,0,0); + + minimizer->SetParameter(0, "vertex0", initialVertex[icand][0], stepSize, initialVertex[icand][0] - limits[0], initialVertex[icand][0] + limits[0]); + minimizer->SetParameter(1, "vertex1", initialVertex[icand][1], stepSize, initialVertex[icand][1] - limits[1], initialVertex[icand][1] + limits[1]); + minimizer->SetParameter(2, "vertex2", initialVertex[icand][2], stepSize, initialVertex[icand][2] - limits[2], initialVertex[icand][2] + limits[2]); + minimizer->SetParameter(3, "vertex3", initialVertex[icand][3], stepSize / fLightSpeed, initialVertex[icand][3] - limits[3] / fLightSpeed, initialVertex[icand][3] + limits[3] / fLightSpeed); + + minimizer->SetParameter(4, "PMTConfiguration", 0, 0, 0, 0); // Useless now + minimizer->SetParameter(5, "nhits", nhits, nhits, nhits - 1, nhits + 1); + minimizer->SetParameter(6, "lowerLimit", lowerLimit, 5e-1, -7, -2); + minimizer->SetParameter(7, "upperLimit", upperLimit, 5e-1, 2, 10); + minimizer->SetParameter(8, "expoSigma", 100, 5, 0, 1e3); + // minimizer->SetParameter(8,"scaling factor",1,1,0,1e3); + // minimizer->SetParameter(8,"signal pe",1,1,0,1e3); + minimizer->SetParameter(9, "directionality", directionality, 1, 0, 2); + minimizer->SetParameter(10, "thread", 0, 0, 0, 0); minimizer->FixParameter(4); minimizer->FixParameter(5); minimizer->FixParameter(6); @@ -1181,383 +736,636 @@ std::vector< std::vector > LEAF::MinimizeVertex(std::vector< std::vector minimizer->FixParameter(8); minimizer->FixParameter(9); minimizer->FixParameter(10); - - minimizer->ExecuteCommand("MIGRAD",Mig,2); - std::vector tCandRecoVtxPos(4,0.); - + minimizer->ExecuteCommand("MIGRAD", Mig, 2); + + std::vector tCandRecoVtxPos(4, 0.); + tCandRecoVtxPos[0] = minimizer->GetParameter(0); tCandRecoVtxPos[1] = minimizer->GetParameter(1); tCandRecoVtxPos[2] = minimizer->GetParameter(2); tCandRecoVtxPos[3] = minimizer->GetParameter(3); - + struct FitPosition hPos; hPos.Vtx = tCandRecoVtxPos; hPos.NLL = 0; - + // Get minuit - minuit=minimizer->GetMinuit(); + minuit = minimizer->GetMinuit(); double fmin, fedm, errdef; - int npar,nparx,istat; - minuit->mnstat(fmin,fedm,errdef,npar,nparx,istat); - + int npar, nparx, istat; + minuit->mnstat(fmin, fedm, errdef, npar, nparx, istat); + hPos.NLL = fmin; - + #ifdef VERBOSE_MIN - if(verbose>=1) std::cout<<"NLL="<= 1) + std::cout << "NLL=" << hPos.NLL << std::endl; #endif tVtxContainer.push_back(hPos); - if(VERBOSE>=2) std::cout<<"Initial vertex position = "<ExecuteCommand("SET PRINTOUT", &p1, 1); // quiet mode minuit->SetErrorDef(1); - //minuit->SetPrintLevel(verbose-1); // quiet mode - if(VERBOSE<2){ - minuit->mnexcm("SET NOWarnings",0,0,err); - } - arglist[0]=2; - minuit->mnexcm("SET STR",arglist,1,err);//set strategy sets the number of derivative estimated. in the present case (2), this is the highest number of it with a counter-balance: very slow! - minimizer->SetFCN(MinuitLikelihood);//here is function to minimize - - double Mig[2]={1e6,1e0};//maxcalls and tolerance - - minimizer->SetParameter(0,"vertex0",0,stepSize ,-limits[0] ,limits[0] ); - minimizer->SetParameter(1,"vertex1",0,stepSize ,-limits[1] ,limits[1] ); - minimizer->SetParameter(2,"vertex2",0,stepSize ,-limits[2] ,limits[2] ); - minimizer->SetParameter(3,"vertex3",0,stepSize/fLightSpeed,-limits[3]/fLightSpeed ,limits[3]/fLightSpeed ); - - minimizer->SetParameter(4,"PMTConfiguration",0,0,0,0); // Useless now - minimizer->SetParameter(5,"nhits",nhits,nhits,nhits-1,nhits+1); - minimizer->SetParameter(6,"lowerLimit",lowerLimit,5e-1,-7,-2); - minimizer->SetParameter(7,"upperLimit",upperLimit,5e-1,2,10); - minimizer->SetParameter(8,"expoSigma",100,5,0,1e3); - minimizer->SetParameter(9,"directionality",directionality,1,0,2); + // minuit->SetPrintLevel(verbose-1); // quiet mode + if (VERBOSE < 2) minuit->mnexcm("SET NOWarnings", 0, 0, err); + arglist[0] = 2; + minuit->mnexcm("SET STR", arglist, 1, err); // set strategy sets the number of derivative estimated. in the present case (2), this is the highest number of it with a counter-balance: very slow! + minimizer->SetFCN(MinuitLikelihood); // here is function to minimize + + double Mig[2] = {1e6, 1e0}; // maxcalls and tolerance + + minimizer->SetParameter(0, "vertex0", 0, stepSize, -limits[0], limits[0]); + minimizer->SetParameter(1, "vertex1", 0, stepSize, -limits[1], limits[1]); + minimizer->SetParameter(2, "vertex2", 0, stepSize, -limits[2], limits[2]); + minimizer->SetParameter(3, "vertex3", 0, stepSize / fLightSpeed, -limits[3] / fLightSpeed, limits[3] / fLightSpeed); + + minimizer->SetParameter(4, "PMTConfiguration", 0, 0, 0, 0); // Useless now + minimizer->SetParameter(5, "nhits", nhits, nhits, nhits - 1, nhits + 1); + minimizer->SetParameter(6, "lowerLimit", lowerLimit, 5e-1, -7, -2); + minimizer->SetParameter(7, "upperLimit", upperLimit, 5e-1, 2, 10); + minimizer->SetParameter(8, "expoSigma", 100, 5, 0, 1e3); + minimizer->SetParameter(9, "directionality", directionality, 1, 0, 2); + minimizer->SetParameter(10, "dir_x", 0, 0, -1, 1); + minimizer->SetParameter(11, "dir_y", 0, 0, -1, 1); + minimizer->SetParameter(12, "dir_z", 0, 0, -1, 1); + minimizer->SetParameter(13, "dir_threshold", 0, 0, 0, 0); minimizer->FixParameter(4); minimizer->FixParameter(5); minimizer->FixParameter(6); minimizer->FixParameter(7); minimizer->FixParameter(8); minimizer->FixParameter(9); - - for(int icand=iStart;icandSetParameter(0,"vertex0",initialVertex[icand][0],stepSize ,initialVertex[icand][0]-limits[0] ,initialVertex[icand][0]+limits[0] ); - minimizer->SetParameter(1,"vertex1",initialVertex[icand][1],stepSize ,initialVertex[icand][1]-limits[1] ,initialVertex[icand][1]+limits[1] ); - minimizer->SetParameter(2,"vertex2",initialVertex[icand][2],stepSize ,initialVertex[icand][2]-limits[2] ,initialVertex[icand][2]+limits[2] ); - minimizer->SetParameter(3,"vertex3",initialVertex[icand][3],stepSize/fLightSpeed ,initialVertex[icand][3]-limits[3]/fLightSpeed ,initialVertex[icand][3]+limits[3]/fLightSpeed ); - - minimizer->SetParameter(4,"PMTConfiguration",0,0,0,0); // Useless now - minimizer->SetParameter(5,"nhits",nhits,nhits,nhits-1,nhits+1); - minimizer->SetParameter(6,"lowerLimit",lowerLimit,5e-1,-7,-2); - minimizer->SetParameter(7,"upperLimit",upperLimit,5e-1,2,10); - minimizer->SetParameter(8,"expoSigma",100,5,0,1e3); - minimizer->SetParameter(9,"directionality",directionality,1,0,2); + for (int icand = iStart; icand < iEnd; icand++) + { + minimizer->SetParameter(0, "vertex0", initialVertex[icand][0], stepSize, initialVertex[icand][0] - limits[0], initialVertex[icand][0] + limits[0]); + minimizer->SetParameter(1, "vertex1", initialVertex[icand][1], stepSize, initialVertex[icand][1] - limits[1], initialVertex[icand][1] + limits[1]); + minimizer->SetParameter(2, "vertex2", initialVertex[icand][2], stepSize, initialVertex[icand][2] - limits[2], initialVertex[icand][2] + limits[2]); + minimizer->SetParameter(3, "vertex3", initialVertex[icand][3], stepSize / fLightSpeed, initialVertex[icand][3] - limits[3] / fLightSpeed, initialVertex[icand][3] + limits[3] / fLightSpeed); + + minimizer->SetParameter(4, "PMTConfiguration", 0, 0, 0, 0); // Useless now + minimizer->SetParameter(5, "nhits", nhits, nhits, nhits - 1, nhits + 1); + minimizer->SetParameter(6, "lowerLimit", lowerLimit, 5e-1, -7, -2); + minimizer->SetParameter(7, "upperLimit", upperLimit, 5e-1, 2, 10); + minimizer->SetParameter(8, "expoSigma", 100, 5, 0, 1e3); + minimizer->SetParameter(9, "directionality", directionality, 1, 0, 2); + + if(fDirection_Filter != nullptr && fDirection_Filter->size() > 2) + { + minimizer->SetParameter(10, "dir_x", (*fDirection_Filter)[0], 5e-1, (*fDirection_Filter)[0] - 1, (*fDirection_Filter)[0] + 1); + minimizer->SetParameter(11, "dir_y", (*fDirection_Filter)[1], 5e-1, (*fDirection_Filter)[1] - 1, (*fDirection_Filter)[1] + 1); + minimizer->SetParameter(12, "dir_z", (*fDirection_Filter)[2], 5e-1, (*fDirection_Filter)[2] - 1, (*fDirection_Filter)[2] + 1); + } + minimizer->SetParameter(13, "dir_threshold", FilterThreshold, 5e-1, FilterThreshold - 1, FilterThreshold + 1); minimizer->FixParameter(4); minimizer->FixParameter(5); minimizer->FixParameter(6); minimizer->FixParameter(7); minimizer->FixParameter(8); minimizer->FixParameter(9); - - minimizer->ExecuteCommand("MIGRAD",Mig,2); + minimizer->FixParameter(10); + minimizer->FixParameter(11); + minimizer->FixParameter(12); + minimizer->FixParameter(13); + + minimizer->ExecuteCommand("MIGRAD", Mig, 2); + + std::vector tCandRecoVtxPos(4, 0.); - std::vector tCandRecoVtxPos(4,0.); - tCandRecoVtxPos[0] = minimizer->GetParameter(0); tCandRecoVtxPos[1] = minimizer->GetParameter(1); tCandRecoVtxPos[2] = minimizer->GetParameter(2); tCandRecoVtxPos[3] = minimizer->GetParameter(3); - + struct FitPosition hPos; hPos.Vtx = tCandRecoVtxPos; hPos.NLL = 0; - + // Get minuit - minuit=minimizer->GetMinuit(); + minuit = minimizer->GetMinuit(); double fmin, fedm, errdef; - int npar,nparx,istat; - minuit->mnstat(fmin,fedm,errdef,npar,nparx,istat); - + int npar, nparx, istat; + minuit->mnstat(fmin, fedm, errdef, npar, nparx, istat); + hPos.NLL = fmin; - + tVtxContainer.push_back(hPos); - if(VERBOSE>=2) std::cout<<"Initial vertex position = "< lock // We don't want simultaneous modification of the common output mtx.lock(); - for ( unsigned int iPos = 0; iPos < tVtxContainer.size(); iPos++ ) { - std::vector vPos(5,0.); + for (unsigned int iPos = 0; iPos < tVtxContainer.size(); iPos++) + { + std::vector vPos(5, 0.); struct FitPosition hPos = tVtxContainer[iPos]; - + vPos[0] = hPos.Vtx[0]; vPos[1] = hPos.Vtx[1]; vPos[2] = hPos.Vtx[2]; vPos[3] = hPos.Vtx[3]; vPos[4] = hPos.NLL; - + fThreadOutput.push_back(vPos); - } delete minimizer; mtx.unlock(); } -struct LEAF::FitterOutput LEAF::MakeFit(const HitCollection* lHitCol, const TimeDelta lTriggerTime, bool bMultiPMT) { +void LEAF::FitVertex(FitterOutput &fOutput, std::vector* fDirection_Filter, float FilterThreshold) +{ + double *tLimits = new double[4]; + int iHitsTotal = fHitCollection->Size(); + std::vector< std::vector > fRecoVtxPosFinal; - fHitCollection = lHitCol; - fTriggerTime = lTriggerTime; - fTimeCorrection = fHitCollection->timestamp - fTriggerTime; + TStopwatch timer; + timer.Reset(); + timer.Start(); - // Get Hits total + std::vector tRecoVtxPosCand = this->SearchVertex_Main(iHitsTotal,fSearchVtxTolerance,false,fSTimePDFLimitsQueueNegative,fSTimePDFLimitsQueuePositive,false); + std::vector> tRecoVtxPos; - //int iHitsTotal = fHitInfo.size(); - int iHitsTotal = fHitCollection->Size(); - - if ( fPositionList.size() == 0 ) { - // Make position lists - this->MakePositionList(); + // std::cout << "time after search in FitVertex is: " << timer.RealTime() << std::endl; + // std::vector SNRList; + + // Extract only the vertex positions from the VtxCandidate structure + for (const auto& cand : tRecoVtxPosCand) + { + tRecoVtxPos.push_back({cand.X, cand.Y, cand.Z, cand.T, cand.NLL}); + // SNRList.push_back(cand.SNR); // Store corresponding SNR } - - bool bDirectionnality = false; - if ( bMultiPMT ) { - bDirectionnality = fUseDirectionality; + + //* was removed for computation time + // if(ContainsTrueVtx(&tRecoVtxPosCand, fTrueVtxPos)) + // fOutput.stepOneContainsTrueVtx = 1; + + if (VERBOSE >= 2) + { + for (unsigned int a = 0; a < tRecoVtxPos.size(); a++) + std::cout << "Candidate vertex [ " << a << " ] = " << tRecoVtxPos[a][0] << " , " << tRecoVtxPos[a][1] << " , " << tRecoVtxPos[a][2] << " , " << tRecoVtxPos[a][3] << ", NLL = " << tRecoVtxPos[a][4] << std::endl; + } + + timer.Stop(); + fOutput.Vtx_Search_ComputeTime = timer.RealTime(); + // std::cout << "SearchVertex took: " << timer.RealTime() << std::endl; + + double dStepSizeFinal = 10; + int iToleranceFinal = 1; + if (fStepByStep) //* Not used and not tested + { + // 2.b. Refined vertex search around candidate found previously + // -> Commented since ages, removed on 2020/11/20 + + // 2.c. More refined vertex search around candidate found previously + double dStepSizeFine2 = 100; + int iToleranceFine2 = 20; + for (int i = 0; i < 4; i++) tLimits[i] = fSearchVtxStep; // particleStart[i]*1e-2;//onversion to meter + std::vector> tRecoVtxPosFine2 = this->SearchVertexFine(tRecoVtxPos, tLimits, dStepSizeFine2, iHitsTotal, fSearchVtxTolerance, iToleranceFine2, VERBOSE, true, false, fSTimePDFLimitsQueueNegative, fSTimePDFLimitsQueuePositive, fUseDirectionality); + ////////////////////////////////////////// + + // 2.d. Final refined vertex search around candidate found previously + double dStepSizeFine3 = 20; + int iToleranceFine3 = 1; + for (int i = 0; i < 4; i++) + tLimits[i] = dStepSizeFine2; // particleStart[i]*1e-2;//onversion to meter + std::vector> tRecoVtxPosFine3 = this->SearchVertexFine(tRecoVtxPosFine2, tLimits, dStepSizeFine3, iHitsTotal, iToleranceFine2, iToleranceFine3, VERBOSE, true, false, fSTimePDFLimitsQueueNegative, fSTimePDFLimitsQueuePositive, fUseDirectionality); + if (VERBOSE == 1) + { + for (int a = 0; a < iToleranceFine3; a++) + std::cout << "Final candidate vertex time = " << tRecoVtxPosFine3[a][0] << ", x = " << tRecoVtxPosFine3[a][1] << ", y = " << tRecoVtxPosFine3[a][2] << ", z=" << tRecoVtxPosFine3[a][3] << std::endl; + } + ////////////////////////////////////////// + + for (int i = 0; i < 4; i++) + tLimits[i] = 100; // dStepSizeFine2;//particleStart[i]*1e-2;//onversion to meter + fRecoVtxPosFinal = this->MinimizeVertex(tRecoVtxPosFine3, tLimits, dStepSizeFinal, iHitsTotal, iToleranceFine3, iToleranceFinal, VERBOSE, true, false, fMinimizeLimitsNegative, fMinimizeLimitsPositive, fUseDirectionality); + if (VERBOSE >= 1) + { + for (int a = 0; a < iToleranceFinal; a++) + { + std::cout << "Final candidate vertex time = " << fRecoVtxPosFinal[a][0] << ", x = " << fRecoVtxPosFinal[a][1] << ", y = " << fRecoVtxPosFinal[a][2] << ", z=" << fRecoVtxPosFinal[a][3] << std::endl; + } + } + } + else + { + timer.Reset(); + timer.Start(); + + for (int i = 0; i < 4; i++) tLimits[i] = 2 * fSearchVtxStep; + + // fRecoVtxPosFinal = this->MinimizeVertex(tRecoVtxPos,tLimits,dStepSizeFinal,iHitsTotal,fSearchVtxTolerance,iToleranceFinal,VERBOSE,true,false,fMinimizeLimitsNegative,fMinimizeLimitsPositive,fUseDirectionality); + fRecoVtxPosFinal = this->MinimizeVertex_Main(tRecoVtxPos, tLimits, dStepSizeFinal, iHitsTotal, fSearchVtxTolerance, iToleranceFinal, VERBOSE, true, false, fMinimizeLimitsNegative, fMinimizeLimitsPositive, fUseDirectionality, fDirection_Filter, FilterThreshold); + + timer.Stop(); + fOutput.Vtx_Minimize_ComputeTime = timer.RealTime(); + // std::cout << "Minimizer took: " << timer.RealTime() << std::endl; } + + fOutput.Vtx[0] = fRecoVtxPosFinal[0][0]; + fOutput.Vtx[1] = fRecoVtxPosFinal[0][1]; + fOutput.Vtx[2] = fRecoVtxPosFinal[0][2]; + fOutput.Vtx[3] = fRecoVtxPosFinal[0][3]; + fOutput.NLL = fRecoVtxPosFinal[0][4]; + // fOutput->NLLR = Likelihoods::GoodnessOfFit(fHitCollection, fRecoVtxPosFinal[0], iHitsTotal, fMinimizeLimitsNegative, fMinimizeLimitsPositive, true, false, fUseDirectionality ); + // fOutput.SNRList = SNRList; + + // fOutput.VtxHitHangles = std::vector(); + + // for (unsigned int ihit = 0; ihit < fHitCollection->Size(); ihit++) + // { + //* compute angle + // Hit lHit = fHitCollection->At(ihit); + // int iPMT = lHit.PMT; + // PMTInfo lPMTInfo = (*fPMTList)[iPMT]; + // std::vector toPMT = std::vector(3); + // for(int j = 0; j < 3; j++) toPMT[j] = lPMTInfo.Position[j] - fRecoVtxPosFinal[0][j]; + // Normalize(toPMT); + // Normalize(fTrueDir); + // double dotP = dot(toPMT, fTrueDir); + // double angle = std::acos(std::max(-1.0, std::min(1.0, dotP))) * 180 / TMath::Pi(); + // fOutput.VtxHitHangles.push_back(angle); + // } +} + +void LEAF::FitEnergy(FitterOutput &fOutput) +{ + int iInTime = 0; + double totalPE = 0; + double totalPECorr = 0; + + TStopwatch timer; + timer.Reset(); + timer.Start(); - //Proceed to the fit. - double * tLimits = new double[4]; + double a = -0.000440808; + double b = 1.71313; + double c = -2.01675; + double d = 1.81755; - struct FitterOutput fOutput; + for (unsigned int ihit = 0; ihit < fHitCollection->Size(); ihit++) + { + // Hit lHit = fHitInfo[ihit]; + Hit lHit = fHitCollection->At(ihit); + + int iPMT = lHit.PMT; + PMTInfo lPMTInfo = (*fPMTList)[iPMT]; + double PMTOrientation[3]; + double dirVector[3]; + for(int j = 0; j<3; j++) PMTOrientation[j] = lPMTInfo.Orientation[j]; + for(int j = 0; j<3; j++) dirVector[j] = -(fTrueVtxPos[j] - lPMTInfo.Position[j]); + GeoTools::Normalize(PMTOrientation); + GeoTools::Normalize(dirVector); + double cosPMThitAngle = abs(dirVector[0] * PMTOrientation[0] + dirVector[1] * PMTOrientation[1] + dirVector[2] * PMTOrientation[2]); + double Ftheta = a + b*cosPMThitAngle + c*(pow(cosPMThitAngle,2)) + d*(pow(cosPMThitAngle,3)); + + // totalPECorr += (lHit.Q/(cosPMThitAngle/**sin(PMThitAngle)*/)) * Ftheta; + // std::cout <<" correction : " << Ftheta / cosPMThitAngle << " for cosPMThitAngle = " << cosPMThitAngle << std::endl; + // totalPECorr += lHit.Q / cosPMThitAngle; + totalPECorr += lHit.Q; + // std::cout << "corr: " << Ftheta / cosPMThitAngle << std::endl; + totalPE += lHit.Q; + double residual = ComputeResidualTime(fOutput.Vtx, fOutput.Vtx[3], fHitCollection->At(ihit)); + if (residual > fSTimePDFLimitsQueueNegative && residual < fSTimePDFLimitsQueuePositive) iInTime++; + } - fOutput.Vtx[0] = 0.; - fOutput.Vtx[1] = 0.; - fOutput.Vtx[2] = 0.; - fOutput.Vtx[3] = 0.; + fOutput.InTime = iInTime; + // double ProportionalFactor = 0.0881523; //* computed experimentally by looking at the steepness of estimated energy vs true energy + // double plus = 4.826; + // ProportionalFactor = 1; + // double energyEstimate = totalPECorr * ProportionalFactor / plus; + double energyEstimate = totalPECorr * 0.111247 - 12.970173; + fOutput.TotalCharge = totalPE; + fOutput.Energy = energyEstimate; + + timer.Stop(); + fOutput.Energy_Fit_ComputeTime = timer.RealTime(); +} + +std::vector LEAF::SearchVertexAndDir(FitterOutput &fOutput) +{ + TStopwatch timer; + timer.Reset(); + timer.Start(); + + std::vector Candidates = std::vector(); + + std::vector tVtxCandidates = this->SearchVertex_Main(fHitCollection->Size(),fSearchVtxTolerance,false,fSTimePDFLimitsQueueNegative,fSTimePDFLimitsQueuePositive,false); + + fOutput.stepOneContainsTrueVtx = (int)ContainsTrueVtx(&tVtxCandidates, fTrueVtxPos); + + for(int i=0; i<(int)tVtxCandidates.size(); i++) + { + std::vector vtx = {tVtxCandidates[i].X, tVtxCandidates[i].Y, tVtxCandidates[i].Z, tVtxCandidates[i].T}; + std::vector tDirCandidates = this->SearchDirection(vtx, DirTolerance); + for(int j=0; j<(int)tDirCandidates.size(); j++) + { + JointFitCandidate candidate; + candidate.VtxPart = tVtxCandidates[i]; + candidate.DirPart = tDirCandidates[j]; + // double DNLL = Likelihoods::Dir_NLL(fHitCollection ,candidate.Vtx, theta, phi, fHitCollection->Size()); + // double VNLL = Likelihoods::Vertex_Score(fHitCollection, candidate.Vtx, fHitCollection->Size(), fMinimizeLimitsNegative, fMinimizeLimitsPositive, true, false, fUseDirectionality); + candidate.NLL = candidate.DirPart.DNLL * candidate.VtxPart.NLL; + Candidates.push_back(candidate); + } + } + + std::sort(Candidates.begin(), Candidates.end(), Likelihoods::CompareByNLL{}); + + while ((int)Candidates.size() > fSearchVtxTolerance) Candidates.pop_back(); + + std::vector dirCandidates = std::vector(); + for(int i=0; i<(int)Candidates.size(); i++) dirCandidates.push_back(Candidates[i].DirPart); + fOutput.stepOneContainsTrueDir = (int)ContainsTrueDir(&dirCandidates, fTrueDir); + + timer.Stop(); + fOutput.Vtx_Search_ComputeTime = timer.RealTime(); + + return Candidates; +} + +JointFitCandidate LEAF::MinimizeVertexAndDir(std::vector initialCandidates, FitterOutput &fOutput) +{ + TStopwatch timer; + timer.Reset(); + timer.Start(); + + double stepSize = 10; + int nhits = fHitCollection->Size(); + double *limits = new double[4]; + for (int i = 0; i < 4; i++) limits[i] = 2 * fSearchVtxStep; + + std::vector tContainer; + + TFitter *minimizer = new TFitter(12); + TMinuit *minuit = minimizer->GetMinuit(); + + double arglist[20]; + int err = 0; + double p1 = VERBOSE - 1; + minimizer->ExecuteCommand("SET PRINTOUT", &p1, 1); // quiet mode + minuit->SetErrorDef(1); + if (VERBOSE < 2) minuit->mnexcm("SET NOWarnings", 0, 0, err); + arglist[0] = 2; + minuit->mnexcm("SET STR", arglist, 1, err); // set strategy sets the number of derivative estimated. in the present case (2), this is the highest number of it with a counter-balance: very slow! + minimizer->SetFCN(MinuitJointNLL); // here is function to minimize + + double Mig[2] = {1e6, 1e0}; // maxcalls and tolerance + + minimizer->SetParameter(0, "vertex0", 0, stepSize, -limits[0], limits[0]); + minimizer->SetParameter(1, "vertex1", 0, stepSize, -limits[1], limits[1]); + minimizer->SetParameter(2, "vertex2", 0, stepSize, -limits[2], limits[2]); + minimizer->SetParameter(3, "vertex3", 0, stepSize / fLightSpeed, -limits[3] / fLightSpeed, limits[3] / fLightSpeed); + minimizer->SetParameter(4, "theta", 0, 2*M_PI/180, 0.0, TMath::Pi()); + minimizer->SetParameter(5, "phi", 0, 2*M_PI/180, -TMath::Pi(), TMath::Pi()); + + minimizer->SetParameter(6, "PMTConfiguration", 0, 0, 0, 0); // Useless now + minimizer->SetParameter(7, "nhits", nhits, nhits, nhits - 1, nhits + 1); + minimizer->SetParameter(8, "lowerLimit", fSTimePDFLimitsQueueNegative, 5e-1, -7, -2); + minimizer->SetParameter(9, "upperLimit", fSTimePDFLimitsQueuePositive, 5e-1, 2, 10); + minimizer->SetParameter(10, "expoSigma", 100, 5, 0, 1e3); + minimizer->SetParameter(11, "thread", 0, 0, 0, 0); + minimizer->FixParameter(6); + minimizer->FixParameter(7); + minimizer->FixParameter(8); + minimizer->FixParameter(9); + minimizer->FixParameter(10); + minimizer->FixParameter(11); + + for (int icand = 0; icand < (int)initialCandidates.size(); icand++) + { + minimizer->SetParameter(0, "vertex0", initialCandidates[icand].VtxPart.X, stepSize, initialCandidates[icand].VtxPart.X - limits[0], initialCandidates[icand].VtxPart.X + limits[0]); + minimizer->SetParameter(1, "vertex1", initialCandidates[icand].VtxPart.Y, stepSize, initialCandidates[icand].VtxPart.Y - limits[1], initialCandidates[icand].VtxPart.Y + limits[1]); + minimizer->SetParameter(2, "vertex2", initialCandidates[icand].VtxPart.Z, stepSize, initialCandidates[icand].VtxPart.Z - limits[2], initialCandidates[icand].VtxPart.Z + limits[2]); + minimizer->SetParameter(3, "vertex3", initialCandidates[icand].VtxPart.T, stepSize / fLightSpeed, initialCandidates[icand].VtxPart.T - limits[3] / fLightSpeed, initialCandidates[icand].VtxPart.T + limits[3] / fLightSpeed); + minimizer->SetParameter(4, "theta", initialCandidates[icand].DirPart.theta, 2*M_PI/180, 0.0, TMath::Pi()); + minimizer->SetParameter(5, "phi", initialCandidates[icand].DirPart.phi, 2*M_PI/180, -TMath::Pi(), TMath::Pi()); + + minimizer->ExecuteCommand("MIGRAD", Mig, 2); + + JointFitCandidate tCandReco; + + tCandReco.VtxPart.X = minimizer->GetParameter(0); + tCandReco.VtxPart.Y = minimizer->GetParameter(1); + tCandReco.VtxPart.Z = minimizer->GetParameter(2); + tCandReco.VtxPart.T = minimizer->GetParameter(3); + tCandReco.DirPart.theta = minimizer->GetParameter(4); + tCandReco.DirPart.phi = minimizer->GetParameter(5); + + // Get minuit + double fmin, fedm, errdef; + int npar, nparx, istat; + minuit->mnstat(fmin, fedm, errdef, npar, nparx, istat); + tCandReco.NLL = fmin; + + tContainer.push_back(tCandReco); + } + + std::sort(tContainer.begin(), tContainer.end(), Likelihoods::CompareByNLL{}); + + delete minimizer; + + timer.Stop(); + fOutput.Vtx_Minimize_ComputeTime = timer.RealTime(); + + return tContainer[0]; +} + +void LEAF::LoadHitCollection(const HitCollection *lHitCol, const TimeDelta lTriggerTime, bool bMultiPMT) +{ + fHitCollection = lHitCol; + fTriggerTime = lTriggerTime; + fTimeCorrection = fHitCollection->timestamp - fTriggerTime; + if (fPositionList.size() == 0) MakePositionList(); + fUseDirectionality = bMultiPMT ? fUseDirectionality : false; +} + +struct FitterOutput LEAF::NewOutuput() +{ + struct FitterOutput fOutput; + // fOutput.SNRList = std::vector(); + fOutput.Vtx = std::vector(4); + fOutput.Dir = std::vector(3); fOutput.NLL = 0.; - + fOutput.DNLL = 0.; + fOutput.Energy = 0.; + fOutput.TotalCharge = 0.; fOutput.InTime = 0; - fOutput.True_NLLDiff = 0; fOutput.True_TimeDiff = 0; fOutput.True_TistDiff = 0; + + // fOutput.VtxHitHangles = std::vector(); + return fOutput; +} + +struct FitterOutput LEAF::MakeSequentialFit(const HitCollection *lHitCol, const TimeDelta lTriggerTime, bool bMultiPMT) +{ + TStopwatch timer; + timer.Reset(); + timer.Start(); + + LoadHitCollection(lHitCol, lTriggerTime, bMultiPMT); - fLastLowerLimit = 0.; - fLastUpperLimit = 0.; - - if ( true ) { + struct FitterOutput fOutput = NewOutuput(); + if(fHitCollection->Size() <= 0) return fOutput; - if(iHitsTotal>0){ - - TStopwatch timer; - timer.Reset(); - timer.Start(); - - //2.a. Coarse vertex search from scratch - double dStepSize=fSearchVtxStep;//in centimeters - //std::vector< std::vector > tRecoVtxPos = this->SearchVertex(iHitsTotal,fSearchVtxTolerance,false,fSTimePDFLimitsQueueNegative,fSTimePDFLimitsQueuePositive,false); - std::vector< std::vector > tRecoVtxPos = this->SearchVertex_Main(iHitsTotal,fSearchVtxTolerance,false,fSTimePDFLimitsQueueNegative,fSTimePDFLimitsQueuePositive,false); - - - if(VERBOSE>=2){ - for(unsigned int a=0;a Commented since ages, removed on 2020/11/20 - - //////////////////////////////////////// - - //2.c. More refined vertex search around candidate found previously - double dStepSizeFine2=100; - int iToleranceFine2=20; - for(int i=0;i<4;i++) tLimits[i] = dStepSize;//particleStart[i]*1e-2;//onversion to meter - std::vector< std::vector > tRecoVtxPosFine2 = this->SearchVertexFine(tRecoVtxPos,tLimits,dStepSizeFine2,iHitsTotal,fSearchVtxTolerance,iToleranceFine2,VERBOSE,true,false,fSTimePDFLimitsQueueNegative,fSTimePDFLimitsQueuePositive,bDirectionnality); - ////////////////////////////////////////// - - //2.d. Final refined vertex search around candidate found previously - double dStepSizeFine3=20; - int iToleranceFine3=1; - for(int i=0;i<4;i++) tLimits[i] = dStepSizeFine2;//particleStart[i]*1e-2;//onversion to meter - std::vector< std::vector > tRecoVtxPosFine3 = this->SearchVertexFine(tRecoVtxPosFine2,tLimits,dStepSizeFine3,iHitsTotal,iToleranceFine2,iToleranceFine3,VERBOSE,true,false,fSTimePDFLimitsQueueNegative,fSTimePDFLimitsQueuePositive,bDirectionnality); - if(VERBOSE == 1){ - for(int a=0; a Commented since ages, removed on 2020/11/20 - - - int iInTime = 0; - - for(int ihit = 0; ihit < iHitsTotal; ihit++){ - //Hit lHit = fHitInfo[ihit]; - Hit lHit = fHitCollection->At(ihit); - - int iPMT = lHit.PMT; - PMTInfo lPMTInfo = (*fPMTList)[iPMT]; - - double hitTime = (fTimeCorrection + lHit.T) / TimeDelta::ns; - double distance = Astro_GetDistance(lPMTInfo.Position,fRecoVtxPosFinal[0]); - double tof = distance / fLightSpeed; - double residual = hitTime - tof - fRecoVtxPosFinal[0][3]; - if(residual > fSTimePDFLimitsQueueNegative && residual < fSTimePDFLimitsQueuePositive){ - iInTime++; - } - } - - // Fill Output - fOutput.Vtx[0] = fRecoVtxPosFinal[0][0]; - fOutput.Vtx[1] = fRecoVtxPosFinal[0][1]; - fOutput.Vtx[2] = fRecoVtxPosFinal[0][2]; - fOutput.Vtx[3] = (fRecoVtxPosFinal[0][3] - fTimeCorrection) / TimeDelta::ns; - fOutput.NLL = fRecoVtxPosFinal[0][4]; - - fOutput.InTime = iInTime; - - //std::cout << " Final vertex " << fOutput.Vtx[0] << " " << fOutput.Vtx[1] << " " << fOutput.Vtx[2] << " Time " << fOutput.Vtx[3] << std::endl; - } - } + //* Vertex + FitVertex(fOutput); + + //* Goodness of fit + fOutput.NLLR = Likelihoods::GoodnessOfFit(fHitCollection, fOutput.Vtx, fHitCollection->Size(), fMinimizeLimitsNegative, fMinimizeLimitsPositive, true, false, fUseDirectionality); //* Goodness of Fit + + //* Direction + FitDirectionQuick(fOutput.Vtx, fOutput); + FitDirection(fOutput.Vtx, fOutput, fHitCollection->Size(), true); + fOutput.MyDir = fOutput.Dir; + fOutput.myDNLL = fOutput.DNLL; + FitDirection(fOutput.Vtx, fOutput, fHitCollection->Size(), false); + // FitDirectionQuick_bis(fOutput.Vtx, fOutput); + + //* Energy + FitEnergy(fOutput); + + timer.Stop(); + fOutput.Leaf_ComputeTime = timer.RealTime(); + + return fOutput; +} + +struct FitterOutput LEAF::MakeJointFit(const HitCollection *lHitCol, const TimeDelta lTriggerTime, bool bMultiPMT) +{ + LoadHitCollection(lHitCol, lTriggerTime, bMultiPMT); + + TStopwatch timer; + timer.Reset(); + timer.Start(); + + struct FitterOutput fOutput = NewOutuput(); + if(fHitCollection->Size() <= 0) return fOutput; + + std::vector Candidates = SearchVertexAndDir(fOutput); + JointFitCandidate finalCandidate = MinimizeVertexAndDir(Candidates, fOutput); + + fOutput.Vtx = {finalCandidate.VtxPart.X, finalCandidate.VtxPart.Y, finalCandidate.VtxPart.Z, finalCandidate.VtxPart.T}; + fOutput.Dir = PolarToCartesianNorm({finalCandidate.DirPart.theta, finalCandidate.DirPart.phi}); + fOutput.NLL = finalCandidate.NLL; + fOutput.DNLL = Likelihoods::Dir_NLL(fHitCollection, fOutput.Vtx, finalCandidate.DirPart.theta, finalCandidate.DirPart.phi , fHitCollection->Size()); + + fOutput.NLLR = Likelihoods::GoodnessOfFit(fHitCollection, fOutput.Vtx, fHitCollection->Size(), fMinimizeLimitsNegative, fMinimizeLimitsPositive, true, false, fUseDirectionality); //* Goodness of Fit + //* Energy + FitEnergy(fOutput); + + timer.Stop(); + fOutput.Leaf_ComputeTime = timer.RealTime(); + return fOutput; } diff --git a/leaf/LEAF.hh b/leaf/LEAF.hh index f9f07ad..566831e 100644 --- a/leaf/LEAF.hh +++ b/leaf/LEAF.hh @@ -30,7 +30,6 @@ // Most of default parameters of method 1 and 2 can be set in LEAF.cc Init() method /*****************************************************************************************************/ - #ifndef LEAF_hh #define LEAF_hh @@ -42,10 +41,10 @@ #include #include -//WCSim Headers +//* WCSim Headers #include "WCSimRootGeom.hh" -//ROOT Headers +//* ROOT Headers #include "TFile.h" #include "TFitter.h" #include "TH1D.h" @@ -59,241 +58,92 @@ #include "TSpline.h" #include "TPaletteAxis.h" -//DataModel informations +//* DataModel informations #include "Geometry.h" #include "HitCollection.h" -// Number of PMT configuration: -#define NPMT_CONFIGURATION 2 - -// Hit definition: -#define NormalPMT 0 // B&L hit -#define MiniPMT 1 // 3" PMT hit -#define AllPMT 2 - -#define VERBOSE 0 - -// Verbose level in functions: -#undef VERBOSE_VTX // In SearchVertex -#undef VERBOSE_NLL - -#define VTX_X 0 -#define VTX_Y 1 -#define VTX_Z 2 -#define VTX_T 3 - -#define N_THREAD 12 // Default for sukap - -// Likelihood: -void MinuitLikelihood(int& nDim, double * gout, double & NLL, double par[], int flg); -void MinimizeVertex_CallThread( int iStart, int iIte, - std::vector< std::vector > initialVertex, double * limits, double stepSize, int nhits, - int nCandidates, int tolerance, int verbose, bool likelihood, bool average, - double lowerLimit, double upperLimit, int directionality); - -void SearchVertex_CallThread( int iStart, int iIte, - int nhits,int tolerance,bool likelihood,double lowerLimit, double upperLimit,int directionality); - -inline bool SortOutputVector ( const std::vector& v1, const std::vector& v2 ) { - return v1[4] < v2[4]; -} - -class LEAF { +//* LEAF Headers +#include "Likelihoods.hh" +#include "LeafInputs.hh" +#include "LeafConfig.hh" +#include "LeafSplines.hh" +#include "LeafUtility.hh" +#include "LeafDefinitions.hh" +class LEAF +{ public: static LEAF* GetME(); void DeleteME(); void Initialize(const Geometry* lGeometry); - void SetNThread(int iThread=N_THREAD) { fThread=iThread; } - - struct FitterOutput { - - double Vtx[4]; - double NLL; - int InTime; - - double True_NLLDiff; - double True_TimeDiff; - double True_TistDiff; - }; - + //* Fitter Main Method. Process the whole fit. + struct FitterOutput MakeSequentialFit(const HitCollection* lHitCol, const TimeDelta lTriggerTime, bool bMultiPMT=true); - // Fitter Main Method. Process the whole fit. - struct FitterOutput MakeFit(const HitCollection* lHitCol, const TimeDelta lTriggerTime, bool bMultiPMT=true); + //* slower + struct FitterOutput MakeJointFit(const HitCollection* lHitCol, const TimeDelta lTriggerTime, bool bMultiPMT=true); + + void SearchVertex_thread(int iStart, int iIte, int nhits, int tolerance=1, bool likelihood=false, double lowerLimit=fSTimePDFLimitsQueueNegative, double upperLimit=fSTimePDFLimitsQueuePositive, int directionality=false); - // NLL - // Calculate likelihood function based on input PDF. For now, the function is based on time residuals. - double FindNLL_Likelihood(std::vector vertexPosition, int nhits, double lowerLimit, double upperLimit, bool killEdges, bool scaleDR, int directionality); - // Calculate likelihood without using PDF, but just using hits within a given timing window (hits "in-time"). - double FindNLL_NoLikelihood(std::vector vertexPosition, int nhits, double lowerLimit, double upperLimit, bool killEdges, bool scaleDR, int directionality); - // Contain the two functions above in one function, where usage of PDF or not can be set through the flag: likelihood = true/false - double FindNLL( std::vector vertexPosition,int nhits, bool likelihood, int verbose, double lowerLimit, double upperLimit, - bool killEdges=false, bool scaleDR=false, int directionality=false); - - void SearchVertex_thread( int iStart, int iIte, - int nhits, - int tolerance=1, bool likelihood=false, - double lowerLimit=fSTimePDFLimitsQueueNegative, double upperLimit=fSTimePDFLimitsQueuePositive, int directionality=false); - - - void MinimizeVertex_thread( int iStart, int iIte, - std::vector< std::vector > initialVertex, double * limits, double stepSize, int nhits, - int nCandidates = 1, int tolerance = 1, int verbose=0, bool likelihood=false, bool average=false, - double lowerLimit=fSTimePDFLimitsQueueNegative, double upperLimit=fSTimePDFLimitsQueuePositive, int directionality = true); + void MinimizeVertex_thread(int iStart, int iIte, std::vector< std::vector > initialVertex, double * limits, double stepSize, int nhits, int nCandidates = 1, int tolerance = 1, int verbose=0, bool likelihood=false, bool average=false, double lowerLimit=fSTimePDFLimitsQueueNegative, double upperLimit=fSTimePDFLimitsQueuePositive, int directionality = true, std::vector* fDirection_Filter = nullptr, float FilterThreshold = 180); private: - + LEAF(); ~LEAF(); - - void Init(); - void LoadSplines(); - void LoadPMTInfo(); - void MakeMPMTReferencial(int iPMT); - void MakePositionList(); - - double GetDistanceOld(std::vector A, std::vector B); - - // Spline functions: - double SplineIntegral (TSpline3 * s, double start, double end, double stepSize=5e-1); - double SplineIntegralAndSubstract (TSpline3 * s0, TSpline3 * s1, double start, double end, double stepSize=5e-1); - double SplineIntegralExpo (TSpline3 * s, double start, double end, double sigma, double stepSize=5e-1); - - // NLL - void MakeEventInfo(double lowerLimit, double upperLimit); - double FindNLLDirectionality(std::vector vertexPosition, int nhits, int verbose, double lowerLimit, double upperLimit); - - - double FindDirectionTheta(std::vector vertex,int tubeNumber, int verbose); + struct FitterOutput NewOutuput(); - // SearchVertex functions: - std::vector< std::vector > SearchVertex( - int nhits, - int tolerance=1, bool likelihood=false, - double lowerLimit=fSTimePDFLimitsQueueNegative, double upperLimit=fSTimePDFLimitsQueuePositive, int directionality=false); - - std::vector< std::vector > SearchVertex_Main( - int nhits, - int tolerance=1, bool likelihood=false, - double lowerLimit=fSTimePDFLimitsQueueNegative, double upperLimit=fSTimePDFLimitsQueuePositive, int directionality=false); - - - std::vector< std::vector > SearchVertexFine( - std::vector< std::vector > initialVertex, double * limits, double stepSize, int nhits, - int nCandidates = 1, int tolerance = 1, int verbose=0, bool likelihood=false, bool average=false, - double lowerLimit=fSTimePDFLimitsQueueNegative, double upperLimit=fSTimePDFLimitsQueuePositive, int directionality=false); - - // Minimize function: - std::vector< std::vector > MinimizeVertex( - std::vector< std::vector > initialVertex, double * limits, double stepSize, int nhits, - int nCandidates = 1, int tolerance = 1, int verbose=0, bool likelihood=false, bool average=false, - double lowerLimit=fSTimePDFLimitsQueueNegative, double upperLimit=fSTimePDFLimitsQueuePositive, int directionality = true); - - std::vector< std::vector > MinimizeVertex_Main( - std::vector< std::vector > initialVertex, double * limits, double stepSize, int nhits, - int nCandidates = 1, int tolerance = 1, int verbose=0, bool likelihood=false, bool average=false, - double lowerLimit=fSTimePDFLimitsQueueNegative, double upperLimit=fSTimePDFLimitsQueuePositive, int directionality = true); + void LoadHitCollection(const HitCollection *lHitCol, const TimeDelta lTriggerTime, bool bMultiPMT); + //* VERTEX + void FitVertex(FitterOutput &fOutput, std::vector* fDirection_Filter = nullptr, float FilterThreshold = 180); - void VectorVertexPMT( std::vector vertex, int iPMT, double* dAngles ); - - // Variables: - private: - - struct FitPosition { - std::vector Vtx; - double NLL; - }; - struct SortingNLL { - bool operator()(FitPosition const &a, FitPosition const &b) const { - return a.NLL < b.NLL; - } - }; - - struct EventInfo { - int hits; - double SignaloverNoise; - double NoiseIntegral; - double SignalIntegral; - }; + //? Coarse GRID Search Raw version, not parallelized, slow + std::vector< std::vector > SearchVertex(int nhits, int tolerance=1, bool likelihood=false, double lowerLimit=fSTimePDFLimitsQueueNegative, double upperLimit=fSTimePDFLimitsQueuePositive, int directionality=false); - static LEAF* myFitter; - - // Spline - TSpline3 * fSplineTimePDFQueue[NPMT_CONFIGURATION]; - TSpline3 * fSplineTimePDFDarkRate[NPMT_CONFIGURATION]; + //? Coarse GRID Search Parallelized version + std::vector SearchVertex_Main(int nhits, int tolerance, bool likelihood, double lowerLimit, double upperLimit, int directionality); - // Histo - TGraph2D * gPMTDirectionality_2D[NPMT_CONFIGURATION][HKAA::kmPMT_Groups]; - - // TF1 - TF1 * fDistResponsePMT[NPMT_CONFIGURATION]; + //? Coarse GRID Search used for step by step vertex fit, written in 2023, not used and not tested + std::vector< std::vector > SearchVertexFine(std::vector< std::vector > initialVertex, double * limits, double stepSize, int nhits, int nCandidates = 1, int tolerance = 1, int verbose=0, bool likelihood=false, bool average=false, double lowerLimit=fSTimePDFLimitsQueueNegative, double upperLimit=fSTimePDFLimitsQueuePositive, int directionality=false); + //? LLH Minimization Raw version, not parallelized, slow + std::vector< std::vector > MinimizeVertex(std::vector< std::vector > initialVertex, double * limits, double stepSize, int nhits, int nCandidates = 1, int tolerance = 1, int verbose=0, bool likelihood=false, bool average=false, double lowerLimit=fSTimePDFLimitsQueueNegative, double upperLimit=fSTimePDFLimitsQueuePositive, int directionality = true); - double fDarkRate_dir_proba[NPMT_CONFIGURATION][HKAA::kmPMT_Groups]; + //? Second Step Parallelized version + std::vector< std::vector > MinimizeVertex_Main(std::vector< std::vector > initialVertex, double * limits, double stepSize, int nhits, int nCandidates = 1, int tolerance = 1, int verbose=0, bool likelihood=false, bool average=false, double lowerLimit=fSTimePDFLimitsQueueNegative, double upperLimit=fSTimePDFLimitsQueuePositive, int directionality = true, std::vector* fDirection_Filter = nullptr, float FilterThreshold = 180); - static double fSTimePDFLimitsQueueNegative; - static double fSTimePDFLimitsQueuePositive; - double fSTimePDFLimitsQueueNegative_fullTimeWindow; - double fSTimePDFLimitsQueuePositive_fullTimeWindow; - double fTimeWindowSizeFull; - - double fMinimizeLimitsNegative; - double fMinimizeLimitsPositive; - - double fHitTimeLimitsNegative; - double fHitTimeLimitsPositive; - - - // Inputs: - const Geometry* fGeometry; - const HitCollection* fHitCollection; - TimeDelta fTriggerTime; - TimeDelta fTimeCorrection; - - // PMT Informations - const std::vector *fPMTList; - double fDarkRate_ns[NPMT_CONFIGURATION]; - - // Detector geometry: - double fTankRadius; - double fTankHeight; - double fTankHalfHeight; - double fLightSpeed; - - // Fitter parameters: - double fIntegrationTimeWindow; - bool fStepByStep; - bool fUseDirectionality; - bool fLimit_mPMT; - int fAveraging; - bool fHighEnergy; - double fSearchVtxStep; - double fSearchVtxTolerance; - - // Input/Output - //std::vector fTrueVtxPos; - //std::vector< std::vector > fTrueVtxPosDouble; - double fPDFNorm_fullTimeWindow; - std::vector< std::vector > fRecoVtxPosFinal; - - - - // Random generator - TRandom3 * fRand; - - double fLastLowerLimit; - double fLastUpperLimit; - struct EventInfo fEventInfo[NPMT_CONFIGURATION]; - - int fThread; - - std::vector< std::vector > fPositionList; + //* DIRECTION + std::vector FitDirection(const std::vector& fixedVertexPosition, FitterOutput &fOutput, int nhits, bool searchPrior = false); + + // Determine direction using the sum of all unit vectors from vertex to PMTs. + std::vector FitDirectionQuick(const std::vector& fixedVertexPosition, FitterOutput &fOutput); + + // ? + std::vector FitDirectionQuick_bis(const std::vector& fixedVertexPosition, FitterOutput &fOutput); + + // Coarse GRID Search + std::vector SearchDirection(const std::vector& fixedVertexPosition, int tolerance); + // LLH Minimization Search + DirectionCandidate MinimizeDirection(const std::vector& fixedVertexPosition, std::vector& candidates, int nhits, double *limits, int verbose); + + + + //##### ENERGY FITTER ##### + + void FitEnergy(FitterOutput &fOutput); + + + //##### JOINT FITTER ##### + + std::vector SearchVertexAndDir(FitterOutput &fOutput); + JointFitCandidate MinimizeVertexAndDir(std::vector initialCandidates, FitterOutput &fOutput); + static LEAF* myFitter; + std::vector< std::vector > fThreadOutput; }; diff --git a/leaf/LeafConfig.cc b/leaf/LeafConfig.cc new file mode 100644 index 0000000..d88c5b3 --- /dev/null +++ b/leaf/LeafConfig.cc @@ -0,0 +1,171 @@ +#include "LeafConfig.hh" + +//* Geometry & Constants +double fTankRadius; +double fTankHeight; +double fTankHalfHeight; +double fLightSpeed; + +//* PDFs +double fSTimePDFLimitsQueueNegative = -3.; +double fSTimePDFLimitsQueuePositive = 4.; +double fSTimePDFLimitsQueueNegative_fullTimeWindow = 0.; +double fSTimePDFLimitsQueuePositive_fullTimeWindow = 0.; +double fPDFNorm_fullTimeWindow = 0.; + +double fDirectionPDF_maxResidual = 15.0; +double fDirectionPDF_minResidual = -5.0; + +//* Overall Fit Parameters +double fTimeWindowSizeFull = 1500.; +double fIntegrationTimeWindow = 50.; +bool fStepByStep = false; +bool fUseDirectionality = false; +bool DoubleNLL = false; +bool fLimit_mPMT = true; +int fAveraging = 20; +bool fHighEnergy = false; +double fSearchVtxStep = 300.; +double fSearchVtxTolerance = 60; +bool bDirectionnality = false; + +//* Step 1 Specific Parameters +double fHitTimeLimitsNegative = -5; +double fHitTimeLimitsPositive = -7; + +//* Step 2 Specific Parameters +double fMinimizeLimitsNegative = -700; +double fMinimizeLimitsPositive = 1000; + +//* Direction Fit Parameters +bool DirTakeAll = false; +double theta_step = 5 * TMath::DegToRad(); +double phi_step = 5 * TMath::DegToRad(); +int DirTolerance = 50; + +void InitConfig(const Geometry *lGeometry) +{ + fTankRadius = lGeometry->detector_radius; + fTankHeight = lGeometry->detector_length; + fTankHalfHeight = fTankHeight / 2.; + + fStepByStep = false; // Step by step mode was not test and is not supported by multithreading + fUseDirectionality = false; + fHighEnergy = false; + DoubleNLL = true; + + // If true, we apply same cuts for timing as B&L PMT when searching vertex. + // Otherwise, we do not apply them and use all the hits. + // The latter is particularly useful when directionality is added, as it can kill DR hits. + fLimit_mPMT = true; + + fSTimePDFLimitsQueueNegative = -3.; //-5;//-3;//-1.5;//-3;//-10; + fSTimePDFLimitsQueuePositive = 4.; // 10;//4;//3;//4;//500; + fSTimePDFLimitsQueueNegative_fullTimeWindow = 0.; //-1.5;//-3;//-10; + fSTimePDFLimitsQueuePositive_fullTimeWindow = 0.; // 3;//4;//500; + fTimeWindowSizeFull = 1500; + fAveraging = 20; // Number of points we used to average the vertex on. + + fMinimizeLimitsNegative = -700; + fMinimizeLimitsPositive = 1000; + + fHitTimeLimitsNegative = -5; //-10;//-5;//-5 //-5 + fHitTimeLimitsPositive = 7; // 15;//7;//7 //7 + + fSearchVtxStep = 300; // in cm, the step size for coarse grid search + + fSearchVtxTolerance = 60; // Number of candidate vertex that are kept after coarse grid search + + fIntegrationTimeWindow = 50; // in ns //? never Used ??? + DirTakeAll = false; + theta_step = 5 * TMath::DegToRad(); + phi_step = 5 * TMath::DegToRad(); + DirTolerance = 50; + + //* Compute light speed: + float fCVacuum = 3e8 * 1e2 / 1e9; // speed of light, in centimeter per ns. + float fNIndex = 1.373; // 1.385;//1.373;//refraction index of water + fLightSpeed = fCVacuum / fNIndex; + + fPDFNorm_fullTimeWindow = 0.; + + std::map envVariables = + { + {"SearchVtxStep", &fSearchVtxStep}, + {"SearchVtxTolerance", &fSearchVtxTolerance}, + {"HitTimeLimitsNegative", &fHitTimeLimitsNegative}, + {"HitTimeLimitsPositive", &fHitTimeLimitsPositive}, + {"MinimizeLimitsNegative", &fMinimizeLimitsNegative}, + {"MinimizeLimitsPositive", &fMinimizeLimitsPositive}, + {"IntegrationTimeWindow", &fIntegrationTimeWindow}, + {"TimeWindowSizeFull", &fTimeWindowSizeFull}, + {"STimePDFLimitsQueueNegative", &fSTimePDFLimitsQueueNegative}, + {"STimePDFLimitsQueuePositive", &fSTimePDFLimitsQueuePositive}, + {"STimePDFLimitsQueueNegative_fullTimeWindow", &fSTimePDFLimitsQueueNegative_fullTimeWindow}, + {"STimePDFLimitsQueuePositive_fullTimeWindow", &fSTimePDFLimitsQueuePositive_fullTimeWindow}, + {"DirectionPDF_maxResidual", &fDirectionPDF_maxResidual}, + {"DirectionPDF_minResidual", &fDirectionPDF_minResidual}, + {"ThetaStep", &theta_step}, + {"PhiStep", &phi_step} + }; + + std::map envIntVariables = + { + {"Averaging", &fAveraging}, + {"DirTolerance", &DirTolerance} + }; + + std::map envBoolVariables = + { + {"Limit_mPMT", &fLimit_mPMT}, + {"UseDirectionality", &fUseDirectionality}, + {"StepByStep", &fStepByStep}, + {"HighEnergy", &fHighEnergy}, + {"DoubleNLL", &DoubleNLL}, + {"DirTakeAll", &DirTakeAll}, + }; + + for (const auto& [envName, variable] : envVariables) + { + const char* envValue = std::getenv(envName.c_str()); + if (envValue != nullptr) + { + *variable = std::stod(envValue); + std::cout << envName << " is set to " << *variable << std::endl; + } + else std::cout << "Environment variable " << envName << " is not set. Using default value : "<< *variable << std::endl; + } + + for (const auto& [envName, variable] : envIntVariables) + { + const char* envValue = std::getenv(envName.c_str()); + if (envValue != nullptr) + { + *variable = std::stoi(envValue); + std::cout << envName << " is set to " << *variable << std::endl; + } + else std::cout << "Environment variable " << envName << " is not set. Using default value : "<< *variable << std::endl; + } + + for (const auto& [envName, variable] : envBoolVariables) + { + const char* envValue = std::getenv(envName.c_str()); + if (envValue != nullptr) + { + std::string valueStr(envValue); + if (valueStr == "true" || valueStr == "1") + { + *variable = true; + std::cout << envName << " is set to true" << std::endl; + } + else if (valueStr == "false" || valueStr == "0") + { + *variable = false; + std::cout << envName << " is set to false" << std::endl; + } + else std::cout << "Environment variable " << envName << " has an invalid value. Using default value : "<< *variable << std::endl; + } + else std::cout << "Environment variable " << envName << " is not set. Using default value : "<< *variable << std::endl; + } + +} \ No newline at end of file diff --git a/leaf/LeafConfig.hh b/leaf/LeafConfig.hh new file mode 100644 index 0000000..ed9fa39 --- /dev/null +++ b/leaf/LeafConfig.hh @@ -0,0 +1,73 @@ +#pragma once + +#define NPMT_CONFIGURATION 1 + +// Hit definition: +#define NormalPMT 0 // B&L hit +#define MiniPMT 1 // 3" PMT hit +#define AllPMT 2 + +#define VERBOSE 0 + +#undef VERBOSE_VTX // In SearchVertex +#undef VERBOSE_NLL // In Likelihood Computation + +#define VTX_X 0 +#define VTX_Y 1 +#define VTX_Z 2 +#define VTX_T 3 + +#define N_THREAD 12 // Default for sukap + +#include "WCSimRootGeom.hh" +#include "Geometry.h" + +// extern constexpr float fCVacuum = 3e8 * 1e2 / 1e9; // speed of light, in centimeter per ns. +// extern constexpr float fNIndex = 1.373; +// extern constexpr double fLightSpeed = fCVacuum / fNIndex; + +//* Geometry & Constants +extern double fTankRadius; +extern double fTankHeight; +extern double fTankHalfHeight; +extern double fLightSpeed; + +//* PDFs +extern double fPDFNorm_fullTimeWindow; +extern double fSTimePDFLimitsQueueNegative; +extern double fSTimePDFLimitsQueuePositive; +extern double fSTimePDFLimitsQueueNegative_fullTimeWindow; +extern double fSTimePDFLimitsQueuePositive_fullTimeWindow; + +extern double fDirectionPDF_maxResidual; +extern double fDirectionPDF_minResidual; + +//* Overall Fit Parameters +extern double fTimeWindowSizeFull; +extern double fIntegrationTimeWindow; +extern bool fStepByStep; +extern bool fUseDirectionality; +extern bool DoubleNLL; +extern bool fLimit_mPMT; +// extern double KELimitsPos; +// extern double KELimitsNeg; +extern int fAveraging; +extern bool fHighEnergy; +extern double fSearchVtxStep; +extern double fSearchVtxTolerance; + +//* Step 1 Specific Parameters +extern double fHitTimeLimitsNegative; +extern double fHitTimeLimitsPositive; + +//* Step 2 Specific Parameters +extern double fMinimizeLimitsNegative; +extern double fMinimizeLimitsPositive; + +//* Direction Fit Parameters +extern bool DirTakeAll; +extern double theta_step; +extern double phi_step; +extern int DirTolerance; + +void InitConfig(const Geometry *lGeometry); \ No newline at end of file diff --git a/leaf/LeafDefinitions.hh b/leaf/LeafDefinitions.hh new file mode 100644 index 0000000..ec7b257 --- /dev/null +++ b/leaf/LeafDefinitions.hh @@ -0,0 +1,96 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +//* WCSim Headers +#include "WCSimRootGeom.hh" + +//* ROOT Headers +#include "TFile.h" +#include "TFitter.h" +#include "TH1D.h" +#include "TH2D.h" +#include "TGraph2D.h" +#include "TF1.h" +#include "TMath.h" +#include "TMinuit.h" +#include "TObject.h" +#include "TRandom3.h" +#include "TSpline.h" +#include "TPaletteAxis.h" + +//* DataModel informations +#include "Geometry.h" +#include "HitCollection.h" + +struct FitterOutput +{ + int stepOneContainsTrueVtx; + int stepOneContainsTrueDir; + std::vector Vtx; + double NLL; + double DNLL; + double myDNLL; + double NLLR; + + int InTime; + + double True_NLLDiff; + double True_TimeDiff; + double True_TistDiff; + + double Energy; + double TotalCharge; + std::vector Dir; + std::vector SNRList; + + std::vector MyDir; + std::vector Quick_Dir; + + //Add informations about computation time in the output: + double Leaf_ComputeTime;//Total leaf computational time + double Vtx_Search_ComputeTime;//Vertex coarse grid search computational time + double Vtx_Minimize_ComputeTime;//Vertex MINUIT minimisation computational time + double Dir_Quick_Search_ComputeTime;//Direction using unit vectors from vector to PMTT - computational time + double Dir_Search_ComputeTime;//Direction coarse grid search computational time + double Dir_Minimize_ComputeTime;//Direction MINUIT minimisation computational time + double Energy_Fit_ComputeTime;//Energy finder computational time +}; + +struct DirectionCandidate +{ + double theta; + double phi; + double DNLL; +}; + +struct VtxCandidate +{ + double X; + double Y; + double Z; + double T; + double NLL; + double SNR; // SNR value +}; + +struct JointFitCandidate +{ + VtxCandidate VtxPart; + DirectionCandidate DirPart; + // std::vector Vtx; + // std::vector Dir; + double NLL; +}; + +struct FitPosition +{ + std::vector Vtx; + double NLL; +}; diff --git a/leaf/LeafInputs.cc b/leaf/LeafInputs.cc new file mode 100644 index 0000000..bcfc9dc --- /dev/null +++ b/leaf/LeafInputs.cc @@ -0,0 +1,83 @@ +#include "LeafInputs.hh" + +const Geometry* fGeometry; + +const HitCollection* fHitCollection; +TimeDelta fTriggerTime; +TimeDelta fTimeCorrection; + +const std::vector *fPMTList; +double fDarkRate_ns[NPMT_CONFIGURATION]; + +int fThread = N_THREAD; + +TRandom3 * fRand; + +std::vector fTrueVtxPos; +std::vector fTrueDir; + +std::vector< std::vector > fPositionList; + +void InitInputs(const Geometry *lGeometry) +{ + fGeometry = lGeometry; + fPMTList = fGeometry->GetPMTList(); + fThread = N_THREAD; + fRand = new TRandom3(); + + fDarkRate_ns[NormalPMT] = fGeometry->pmt_dark_rate[HKAA::kIDPMT_BnL]; + fDarkRate_ns[MiniPMT] = fGeometry->pmt_dark_rate[HKAA::kIDPMT_3inch]; + MakePositionList(); +} + +void MakePositionList() +{ + // Make position list for a given step size + double dStep = fSearchVtxStep; + fPositionList = std::vector< std::vector >(); + fPositionList.clear(); + + for (double time = -50; time < 50; time += (dStep / fLightSpeed)) + { + for (double radius = 0; radius <= fTankRadius; radius += dStep) + { + double perimeter = 2 * TMath::Pi() * radius; + int numberOfStepsOnCircle = floor(perimeter / dStep); + if (numberOfStepsOnCircle == 0) numberOfStepsOnCircle = 1; + double angleStepOnCircle = 2 * TMath::Pi() / numberOfStepsOnCircle; + + for (double angle = 0; angle <= 2 * TMath::Pi(); angle += angleStepOnCircle) + { + for (double height = -fTankHalfHeight; height <= fTankHalfHeight; height += dStep) + { + std::vector tVtxPosition(4, 0.); // In centimeters + tVtxPosition[0] = radius * TMath::Cos(angle); + tVtxPosition[1] = radius * TMath::Sin(angle); + tVtxPosition[2] = height; + tVtxPosition[3] = time; + fPositionList.push_back(tVtxPosition); + } + } + } + } +} + +void SetTrueVertexInfo(std::vector vtx, double time) +{ + fTrueVtxPos = std::vector(5, 0.); + fTrueVtxPos[0] = vtx[0]; + fTrueVtxPos[1] = vtx[1]; + fTrueVtxPos[2] = vtx[2]; + fTrueVtxPos[3] = time; + fTrueVtxPos[4] = 0.; +} + +void SetTrueDirInfo(std::vector trueDir) +{ + fTrueDir = std::vector(3, 0.); + fTrueDir[0] = trueDir[0]; + fTrueDir[1] = trueDir[1]; + fTrueDir[2] = trueDir[2]; +} + +void SetNThread(int iThread){ fThread=iThread; } \ No newline at end of file diff --git a/leaf/LeafInputs.hh b/leaf/LeafInputs.hh new file mode 100644 index 0000000..a9643fd --- /dev/null +++ b/leaf/LeafInputs.hh @@ -0,0 +1,47 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +//* Root Headers +#include "TRandom3.h" + +//* WCSim Headers +#include "WCSimRootGeom.hh" +#include "Geometry.h" +#include "HitCollection.h" +#include "LeafConfig.hh" + +extern const Geometry* fGeometry; +extern const HitCollection* fHitCollection; +extern TimeDelta fTriggerTime; +extern TimeDelta fTimeCorrection; + +extern const std::vector *fPMTList; +extern double fDarkRate_ns[NPMT_CONFIGURATION]; + +extern std::vector< std::vector > fPositionList; + +extern int fThread; + +extern TRandom3 * fRand; + +extern std::vector fTrueVtxPos; +extern std::vector fTrueDir; + +void InitInputs(const Geometry *lGeometry); + +void MakePositionList(); + +void SetNThread(int iThread=N_THREAD); + +//* True Vertex isn't used for the fit, just to get feedback on how well the fit is at different steps +void SetTrueVertexInfo(std::vector vtx, double time); +void SetTrueDirInfo(std::vector trueDir); + + diff --git a/leaf/LeafSplines.cc b/leaf/LeafSplines.cc new file mode 100644 index 0000000..1abb100 --- /dev/null +++ b/leaf/LeafSplines.cc @@ -0,0 +1,208 @@ +#include "LeafSplines.hh" + +TSpline3 * fSplineTimePDFQueue[NPMT_CONFIGURATION]; +TSpline3 * fSplineTimePDFDarkRate[NPMT_CONFIGURATION]; +TSpline3* fDirectionPDF; +TSpline3* fHitAnglePDF; +TGraph2D * gPMTDirectionality_2D[NPMT_CONFIGURATION][HKAA::kmPMT_Groups]; +TF1 * fDistResponsePMT[NPMT_CONFIGURATION]; + +double fDarkRate_dir_proba[NPMT_CONFIGURATION][HKAA::kmPMT_Groups]; +double fLastLowerLimit; +double fLastUpperLimit; + +double SplineIntegral(TSpline3 *s, double start, double end, double stepSize) +{ + double integral = 0; + for (double i = start; i < end; i += stepSize) + { + integral += stepSize * s->Eval(i + stepSize / 2); + } + return integral; +} +double SplineIntegralAndSubstract(TSpline3 *s0, TSpline3 *s1, double start, double end, double stepSize) +{ + double integral = 0; + for (double i = start; i < end; i += stepSize) + { + integral += stepSize * (s0->Eval(i + stepSize / 2) - s1->Eval(i + stepSize / 2)); + } + return integral; +} + +double SplineIntegralExpo(TSpline3 *s, double start, double end, double sigma, double stepSize) +{ + double integral = 0; + for (double i = start; i < end; i += stepSize) + { + integral += stepSize * s->Eval(i + stepSize / 2) * TMath::Gaus(i + stepSize / 2, 0, sigma); + } + return integral; +} + +void LoadSplines() +{ + std::cout << "Loading splines..." << std::endl; + TFile *fSplines, *fSplines2, *fSplines1, *fSplineAngleFile; + if (fHighEnergy) + { + fSplines = new TFile("${LEAFDIR}/inputs/timePDF_HE.root", "read"); // To generate with my code ProduceWSPlots.c + fSplines2 = fSplines; + } + else + { + std::cout << "spline 1" << std::endl; + fSplines1 = new TFile("${LEAFDIR}/inputs/timePDFNoDR_50000_e10MeV_Hit_720.root","read");//To generate with my code ProduceWSPlots.c + // fSplines = new TFile("${LEAFDIR}/inputs/timePDF_DRnew_Large.root", "read"); // To generate with my code ProduceWSPlots.c + std::cout << "spline 2" << std::endl; + fSplines = new TFile("${LEAFDIR}/inputs/timePDFDR_5000_e10MeV_Hit_fiducial.root","read"); //*TAHA + // fSplines = new TFile("${LEAFDIR}/inputs/timePDF_3M_10T_500000.root","read"); //*NICOLAS + // fSplines = new TFile("${LEAFDIR}/inputs/timePDF_3M_10T.root","read"); + std::cout << "spline 3" << std::endl; + fSplines2 = new TFile("${LEAFDIR}/inputs/timePDF_Directionality_DRnew.root", "read"); // To generate with my code ProduceWSPlots.c + // fSplines2 = new TFile("${LEAFDIR}/inputs/timePDF_Directionality_DRnew.root","read");//To generate with my code ProduceWSPlots.c + std::cout<< "angle spline" << std::endl; + fSplineAngleFile = new TFile("${LEAFDIR}/inputs/Hit_Angle_PDF.root", "read"); + } + + std::cout << "spline file read" << std::endl; + + // Prevent TGraph2D to be append to TFile (this is needed as we are doing multiple copy of TGraph2D) + // For a strange reason sometimes the TGraph2D is see as a TH1 and stored in an open TFile + TH1::AddDirectory(false); + + //Check if the files are opened + if (!fSplines || !fSplines->IsOpen()) + { + std::cerr << "Error: Splines file is not opened." << std::endl; + exit(1); + } + if (!fSplines2 || !fSplines2->IsOpen()) + { + std::cerr << "Error: Splines file is not opened." << std::endl; + exit(1); + } + if (!fSplines1 || !fSplines1->IsOpen()) + { + std::cerr << "Error: Splines file is not opened." << std::endl; + exit(1); + } + if (!fSplineAngleFile || !fSplineAngleFile->IsOpen()) + { + std::cerr << "Error: Angle spline file is not opened." << std::endl; + exit(1); + } + + // std::cout << "splines files loaded" << std::endl; + + int configs = NPMT_CONFIGURATION; + + for(int pmtType=0;pmtTypeGet(Form("splineExpoQueue%d_%d",0,pmtType)); + fSplineTimePDFDarkRate[pmtType] = (TSpline3*) fSplines->Get(Form("splineDR%d_%d",0,pmtType)); + fDirectionPDF = (TSpline3*) fSplines1->Get(Form("ChargeProfileCos_PDF_spline_pmtType0")); + // fDirectionPDF = (TSpline3*) fSplineAngleFile->Get("Hit_Angle_Spline"); + if (!fDirectionPDF) + { + std::cerr << "Error: fDirectionPDF is not loaded properly." << std::endl; + exit(1); + } + fHitAnglePDF = (TSpline3*) fSplineAngleFile->Get("Hit_Angle_Spline"); + // fHitAnglePDF = (TSpline3*) fSplines1->Get(Form("ChargeProfileCos_PDF_spline_pmtType0")); + + std::cout << "process spline" << std::endl; + // std::cout << "ok ?" << std::endl; + fSTimePDFLimitsQueueNegative_fullTimeWindow = fSplineTimePDFQueue[pmtType]->GetXmin(); + fSTimePDFLimitsQueuePositive_fullTimeWindow = fSplineTimePDFQueue[pmtType]->GetXmax(); + //std::cout<<"Min="<IsOpen()) { + std::cerr << "Error: fSplines2 is not opened or is null." << std::endl; + exit(1); + } + + int grpNB = 0; + if (HKAA::kmPMT_Groups > 0) { + grpNB = HKAA::kmPMT_Groups; + } else { + std::cerr << "Error: HKAA::kmPMT_Groups is not properly defined or initialized." << std::endl; + exit(1); + } + + std::cout << "got grnNB" << std::endl; + + for (int pmtGroup = 0; pmtGroup < grpNB; pmtGroup++) { + gPMTDirectionality_2D[pmtType][pmtGroup] = (TGraph2D*) fSplines2->Get(Form("gPMTDirectionality_2D_%d_%d_%d", 0, pmtType, pmtGroup)); + if (!gPMTDirectionality_2D[pmtType][pmtGroup]) { + std::cerr << "Error: Failed to load gPMTDirectionality_2D for pmtType " << pmtType << " and pmtGroup " << pmtGroup << std::endl; + exit(1); + } + } + + std::cout << "dist response" << std::endl; + + fDistResponsePMT[pmtType] = (TF1*) fSplines2->Get(Form("fDistResponsePMT_pmtType%d", pmtType)); + if (!fDistResponsePMT[pmtType]) + { + std::cerr << "Error: Failed to load fDistResponsePMT for pmtType " << pmtType << std::endl; + exit(1); + } + } + std::cout << "Splines loaded." << std::endl; +} + +EventInfo MakeEventInfo(double lowerLimit, double upperLimit, int pmtType) +{ + struct EventInfo fEventInfo; + + if (fLastLowerLimit == lowerLimit && fLastUpperLimit == upperLimit) return fEventInfo; + fLastLowerLimit = lowerLimit; + fLastUpperLimit = upperLimit; + + fEventInfo.hits = 0; + fEventInfo.SignaloverNoise = 0.; + fEventInfo.NoiseIntegral = 0.; + fEventInfo.SignalIntegral = 0.; + + int iHitTotal = fHitCollection->Size(); + + for (int iHit = 0; iHit < iHitTotal; iHit++) + { + // Hit lHit = fHitInfo[iHit]; + Hit lHit = fHitCollection->At(iHit); + + int iPMT = lHit.PMT; + int iType = Astro_GetPMTType(iPMT); + if(pmtType == iType) fEventInfo.hits += 1; + } + + double signalDR, signalPE; + signalDR = fTimeWindowSizeFull * fDarkRate_ns[pmtType]; // Over the whole time window + signalPE = std::max(fEventInfo.hits - signalDR, 0.); // Over the whole time window. + + double signalPETime = signalPE; // I assume that all signal is here. + double signalDRTime = (upperLimit - lowerLimit) * signalDR / fTimeWindowSizeFull; + + fEventInfo.SignaloverNoise = signalPETime / signalDRTime; // Over the + fEventInfo.NoiseIntegral = SplineIntegral(fSplineTimePDFDarkRate[pmtType], lowerLimit, upperLimit); + fEventInfo.SignalIntegral = SplineIntegralAndSubstract(fSplineTimePDFQueue[pmtType], fSplineTimePDFDarkRate[pmtType], lowerLimit, upperLimit); + + if (VERBOSE >= 3) std::cout + << "nhits=" << fEventInfo.hits + << ", DR average=" << signalDR + << ", signal over noise=" << fEventInfo.SignaloverNoise + << ", signal integral=" << fEventInfo.SignalIntegral + << ", DR integral=" << fEventInfo.NoiseIntegral + << ",in integral=" << fEventInfo.SignalIntegral / fEventInfo.NoiseIntegral << std::endl; + + return fEventInfo; +} \ No newline at end of file diff --git a/leaf/LeafSplines.hh b/leaf/LeafSplines.hh new file mode 100644 index 0000000..6d2b783 --- /dev/null +++ b/leaf/LeafSplines.hh @@ -0,0 +1,64 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include "TFile.h" +#include "TFitter.h" +#include "TH1D.h" +#include "TH2D.h" +#include "TGraph2D.h" +#include "TF1.h" +#include "TMath.h" +#include "TMinuit.h" +#include "TObject.h" +#include "TRandom3.h" +#include "TSpline.h" +#include "TPaletteAxis.h" + +//WCSim Headers +#include "WCSimRootGeom.hh" +#include "Geometry.h" +#include "HitCollection.h" +#include "LeafInputs.hh" + + +struct EventInfo +{ + int hits; + double SignaloverNoise; + double NoiseIntegral; + double SignalIntegral; +}; + +// extern bool fHighEnergySplines; + +extern TSpline3 * fSplineTimePDFQueue[NPMT_CONFIGURATION]; +extern TSpline3 * fSplineTimePDFDarkRate[NPMT_CONFIGURATION]; +extern TSpline3* fDirectionPDF; +extern TSpline3* fHitAnglePDF; + +// Histo +extern TGraph2D * gPMTDirectionality_2D[NPMT_CONFIGURATION][HKAA::kmPMT_Groups]; + +// TF1 +extern TF1 * fDistResponsePMT[NPMT_CONFIGURATION]; + +extern double fDarkRate_dir_proba[NPMT_CONFIGURATION][HKAA::kmPMT_Groups]; + +extern double fLastLowerLimit; +extern double fLastUpperLimit; +// double fTimeWindowSizeFull; + +void LoadSplines(); + +double SplineIntegral (TSpline3 * s, double start, double end, double stepSize=5e-1); +double SplineIntegralAndSubstract (TSpline3 * s0, TSpline3 * s1, double start, double end, double stepSize=5e-1); +double SplineIntegralExpo (TSpline3 * s, double start, double end, double sigma, double stepSize=5e-1); + +EventInfo MakeEventInfo(double lowerLimit, double upperLimit, int pmtType); \ No newline at end of file diff --git a/leaf/LeafUtility.cc b/leaf/LeafUtility.cc new file mode 100644 index 0000000..cc2be55 --- /dev/null +++ b/leaf/LeafUtility.cc @@ -0,0 +1,351 @@ +#include "LeafUtility.hh" + +bool ContainsTrueVtx(std::vector* candidates, std::vector fTrueVtxPos) +{ + for (unsigned int i = 0; i < candidates->size(); i++) + { + std::vector vtx = {(*candidates)[i].X, (*candidates)[i].Y, (*candidates)[i].Z, (*candidates)[i].T}; + double distance = Distance3D(vtx, fTrueVtxPos); + if (distance <= fSearchVtxStep) //? should we take step / 2 ? + { + return true; + } + } + return false; +} + +bool ContainsTrueDir(std::vector* candidates, std::vector fTrueDir) +{ + for (unsigned int i = 0; i < candidates->size(); i++) + { + std::vector truePolarDir = CartesianToPolarNorm(fTrueDir); + if(abs((*candidates)[i].theta - truePolarDir[0]) < theta_step && abs((*candidates)[i].phi - truePolarDir[1]) < phi_step) + { + return true; + } + } + return false; +} + + +//* to know if the predicted vertex and the true vertex belong to the same candidate +//? doesn't work because the radius overlap between candidates +bool CanidatesMatch(std::vector>* candidates, std::vector point, std::vector fTrueVtxPos) +{ + for(unsigned int i = 0; i < candidates->size(); i++) + { + if (Distance3D((*candidates)[i], fTrueVtxPos) <= fSearchVtxStep && Distance3D((*candidates)[i], point) <= fSearchVtxStep) return true; + } + return false; +} + +double Distance3D(std::vector point1, std::vector point2) +{ + return sqrt(pow(point1[0] - point2[0], 2) + + pow(point1[1] - point2[1], 2) + + pow(point1[2] - point2[2], 2)); +} + +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; + } +} + +std::vector& Normalize(std::vector& vector) +{ + if (vector.size() < 3) throw std::invalid_argument("Vector must be of size at least 3 to be normalized"); + double length = sqrt(vector[0] * vector[0] + vector[1] * vector[1] + vector[2] * vector[2]); + if (length == 0) throw std::invalid_argument("Cannot normalize a zero-length vector."); + for (int i = 0; i < 3; ++i) vector[i] /= length; + return vector; +} + +double calculateDistance(const std::vector& A, const std::vector& B) +{ + return sqrt(pow(A[0] - B[0], 2) + pow(A[1] - B[1], 2) + pow(A[2] - B[2], 2)); +} + +double dot(const std::vector& A, const std::vector& B) +{ + if (A.size() != B.size()) throw std::invalid_argument("Vectors must be of the same size for dot product."); + double result = 0; + for (size_t i = 0; i < A.size(); ++i) result += A[i] * B[i]; + return result; +} + +std::vector PolarToCartesianNorm(std::vector polarVector) +{ + double theta = polarVector[0]; + double phi = polarVector[1]; + std::vector cartVector = std::vector(3); + cartVector[0] = sin(theta) * cos(phi); + cartVector[1] = sin(theta) * sin(phi); + cartVector[2] = cos(theta); + Normalize(cartVector); + return cartVector; +} + +std::vector CartesianToPolarNorm(const std::vector& cartVector) +{ + std::vector polarVector(2); + double x = cartVector[0]; + double y = cartVector[1]; + double z = cartVector[2]; + + polarVector[0] = acos(z); // theta + polarVector[1] = atan2(y, x); // phi + + return polarVector; +} + +double EstimateBandeWidth(const std::vector& residuals) +{ + double sum = std::accumulate(residuals.begin(), residuals.end(), 0.0); + double mean = sum / residuals.size(); + + double sq_sum = std::inner_product(residuals.begin(), residuals.end(), residuals.begin(), 0.0); + double variance = sq_sum / residuals.size() - mean * mean; + double stddev = std::sqrt(variance); + return 1.06 * stddev * std::pow(residuals.size(), -1.0 / 5.0); +} + +double GaussianKernel(double x, double bandwidth) +{ + return std::exp(-0.5 * x * x / (bandwidth * bandwidth)) / (bandwidth * std::sqrt(2.0 * M_PI)); +} + +double KDE_Estimate(const std::vector& residuals, double t_i, double bandwidth) +{ + double sum = 0.0; + for (const double& t_j : residuals) sum += GaussianKernel(t_i - t_j, bandwidth); + return sum / residuals.size(); +} + +double ComputeResidualTime(std::vector vertexPos, double originTime, Hit lHit) +{ + int iPMT = lHit.PMT; + double HitT = (fTimeCorrection + lHit.T) / TimeDelta::ns; + PMTInfo lPMTInfo = (*fPMTList)[iPMT]; + double distance = Astro_GetDistance(lPMTInfo.Position, vertexPos); + + double tof = distance / fLightSpeed; + return HitT - tof - originTime; +} + +VtxCandidate CreateCandidate(const std::vector& vertex, double lowerLimit, double upperLimit, bool computeSNR) +{ + VtxCandidate out; + + // Copy the vertex info + out.X = vertex[0]; + out.Y = vertex[1]; + out.Z = vertex[2]; + out.T = vertex[3]; + out.NLL = 0.0; + + if(computeSNR) + { + // We'll do a simple approach: count how many hits are in [lower,upper] + // for each PMT type, then subtract the dark rate. + + int nHitsTotal = fHitCollection->Size(); + int inTimeBnL = 0; + int inTimemPMT = 0; + + for(int iHit = 0; iHit < nHitsTotal; iHit++) { + const Hit& hit = fHitCollection->At(iHit); + int iPMT = hit.PMT; + int pmtType = Astro_GetPMTType(iPMT); // 0=BnL, 1=mPMT + + double hitTime = (fTimeCorrection + hit.T) / TimeDelta::ns; + PMTInfo pInfo = (*fPMTList)[iPMT]; + double distance = Astro_GetDistance(pInfo.Position, vertex); + double tof = distance / fLightSpeed; + double residual = hitTime - tof - vertex[3]; + //std::cout << "residual: " << residual << ", " << "lowerlimit and upperlimit " << lowerLimit << ", " << upperLimit << std::endl; + if(residual >= lowerLimit && residual <= upperLimit) + { + if(pmtType == 0) { + // BnL + inTimeBnL++; + } else { + // mPMT + inTimemPMT++; + } + } + } + + double dt = (upperLimit - lowerLimit); + + // BnL: + double darkRateBnL = fTimeWindowSizeFull * fDarkRate_ns[NormalPMT]; + // average dark over dt: + double darkInWindowBnL = (dt / double(fTimeWindowSizeFull)) * darkRateBnL; + double signalBnL = std::max(double(inTimeBnL) - darkInWindowBnL, 0.0); + //std::cout << "inTimeBnL - darkInWindowBnL : " << inTimeBnL << "-" << darkInWindowBnL << std::endl; + out.SNR = (darkInWindowBnL > 1e-9) ? (signalBnL / darkInWindowBnL) : 999.0; + } + return out; +} + +void VectorVertexPMT(std::vector vertex, int iPMT, double *dAngles) +{ + // Guillaume 2020/05/20: + // mPMT referencial is computed once for all PMT in LoadPMTInfo() and MakeMPMTReferencial(iPMT) + // Guillaume 2020/11/20: + // Reference is moved to outside LEAF (HKManager or hk-AstroAnalysis) + + PMTInfo lPMTInfo = (*fPMTList)[iPMT]; + int iPMTTop = lPMTInfo.mPMT_RefTube; + PMTInfo lPMTInfoTop = (*fPMTList)[iPMTTop]; + + // 5. Now we have our referential, we should just calculate the angles of the PMT to vertex position vector in this referential. + // a. calculate the PMT to vertex position vector. + + double dVtx_PMTRef[3]; + + dVtx_PMTRef[0] = vertex[0] - lPMTInfoTop.Position[0]; + dVtx_PMTRef[1] = vertex[1] - lPMTInfoTop.Position[1]; + dVtx_PMTRef[2] = vertex[2] - lPMTInfoTop.Position[2]; + + double dLengthVtx = Astro_GetLength(dVtx_PMTRef); + GeoTools::Normalize(dVtx_PMTRef); + + if (VERBOSE >= 3) + { + std::cout << "Vertex position = " << vertex[0] << ", " << vertex[1] << ", " << vertex[2] << std::endl; + std::cout << "Vertex position from PMT = " << dVtx_PMTRef[0] << ", " << dVtx_PMTRef[1] << ", " << dVtx_PMTRef[2] << std::endl; + } + + // b. Then extract Theta and Phi: + double dCosTheta = Astro_GetScalarProd(dVtx_PMTRef, lPMTInfo.mPMT_RefZ); + double dTheta = TMath::ACos(dCosTheta); + + double dPhi = 0.; + + if ((*fPMTList)[iPMT].mPMT_TubeNum == HKAA::kmPMT_TopID) + { + // Phi is not defined in that case.. + } + else + { + // We know x=cosPhi x sinTheta and y=sinPhi x sinTheta + double dX = Astro_GetScalarProd(dVtx_PMTRef, lPMTInfo.mPMT_RefX); + double dY = Astro_GetScalarProd(dVtx_PMTRef, lPMTInfo.mPMT_RefY); + double dTanPhi = dY / dX; + dPhi = TMath::ATan(dTanPhi); + + // tan is symetric from -pi/2 to +pi/2 + if (dX == 0) + { + if (dY < 0) + dPhi = -TMath::Pi() / 2.; + else + dPhi = TMath::Pi() / 2.; + } + if (dX < 0) + dPhi += TMath::Pi(); // With this, angle become defines between -pi/2 to -pi/2 +2pi. + + // We wish to bring this from 0 to 2pi: + if (dPhi < 0) + dPhi += 2 * TMath::Pi(); + + // Actually, we have a symmetry in Phi between [0,pi] and [pi,2pi]. So, we will just define Phi in [0,pi] modulo pi + if (dPhi > TMath::Pi()) + dPhi = TMath::Pi() - (dPhi - TMath::Pi()); + } + + dAngles[0] = dPhi * 180. / TMath::Pi(); + dAngles[1] = dTheta * 180. / TMath::Pi(); + dAngles[2] = dLengthVtx; + dAngles[3] = 1; // calculateWeight(dLengthVtx,PMTradius[pmtType],Theta,verbose); + + if (VERBOSE >= 3) + { + std::cout << "Angles Phi = " << dAngles[0] << ", Theta = " << dAngles[1] << std::endl; + } +} + +// std::vector ProjectPointToCylinder(std::vector point, double R, double H) +// { +// std::vector dummy; +// return dummy; +// } + +std::vector ProjectPointToCylinder(std::vector point, double R, double H) +{ + double x = point[0]; + double y = point[1]; + double z = point[2]; + double halfHeight = H / 2.0; + + // Distance from the z-axis + double r_xy = std::sqrt(x * x + y * y); + + // --- Project to side surface --- + double x_side, y_side; + if (r_xy == 0.0) + { + x_side = R; + y_side = 0.0; // Arbitrary direction + } + else + { + x_side = x * R / r_xy; + y_side = y * R / r_xy; + } + double z_side = std::clamp(z, -halfHeight, halfHeight); + + // --- Project to top cap (z = +halfHeight) --- + double z_top = halfHeight; + double x_top = x, y_top = y; + double r_top = std::sqrt(x * x + y * y); + if (r_top > R && r_top != 0.0) + { + x_top = x * R / r_top; + y_top = y * R / r_top; + } + + // --- Project to bottom cap (z = -halfHeight) --- + double z_bottom = -halfHeight; + double x_bottom = x, y_bottom = y; + double r_bottom = std::sqrt(x * x + y * y); + if (r_bottom > R && r_bottom != 0.0) + { + x_bottom = x * R / r_bottom; + y_bottom = y * R / r_bottom; + } + + // --- Compute squared distances to each projection --- + auto dist2 = [](double dx, double dy, double dz) + { + return dx * dx + dy * dy + dz * dz; + }; + + double d2_side = dist2(x - x_side, y - y_side, z - z_side); + double d2_top = dist2(x - x_top, y - y_top, z - z_top); + double d2_bottom = dist2(x - x_bottom, y - y_bottom, z - z_bottom); + + // --- Return closest projection --- + if (d2_side <= d2_top && d2_side <= d2_bottom) return {x_side, y_side, z_side}; + else if (d2_top <= d2_bottom) return {x_top, y_top, z_top}; + else return {x_bottom, y_bottom, z_bottom}; +} + +std::vector VectorToHitNorm(const std::vector& from, Hit lHitt) +{ + std::vector toHit(3, 0.0); + std::vector hitPos(3, 0.0); + int iPMT = lHitt.PMT; + PMTInfo lPMTInfo = (*fPMTList)[iPMT]; + for (int j = 0; j < 3; j++) hitPos[j] = lPMTInfo.Position[j]; + for (int j = 0; j < 3; j++) toHit[j] = hitPos[j] - from[j]; + + Normalize(toHit); + + return toHit; +} \ No newline at end of file diff --git a/leaf/LeafUtility.hh b/leaf/LeafUtility.hh new file mode 100644 index 0000000..444de27 --- /dev/null +++ b/leaf/LeafUtility.hh @@ -0,0 +1,59 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +//* DataModel informations +#include "Geometry.h" +#include "HitCollection.h" +#include "LeafConfig.hh" +#include "LeafInputs.hh" +#include "LeafDefinitions.hh" + +inline bool SortOutputVector ( const std::vector& v1, const std::vector& v2 ) +{ + return v1[4] < v2[4]; +} + +double Distance3D(std::vector point1, std::vector point2); + +void Normalize(double a[3]); + +std::vector& Normalize(std::vector& vector); + +double dot(const std::vector& A, const std::vector& B); + +double calculateDistance(const std::vector& A, const std::vector& B); + +std::vector VectorToHitNorm(const std::vector& from, Hit lHitt); + +std::vector PolarToCartesianNorm(std::vector polarVector); + +std::vector CartesianToPolarNorm(const std::vector& cartVector); + +bool CanidatesMatch(std::vector>* candidates, std::vector point, std::vector fTrueVtxPos); + +bool ContainsTrueVtx(std::vector* candidates, std::vector fTrueVtxPos); + +bool ContainsTrueDir(std::vector* candidates, std::vector fTrueDir); + +//* Estimate P_data(t_i) using KDE +double KDE_Estimate(const std::vector& residuals, double t_i, double bandwidth); + +double EstimateBandeWidth(const std::vector& residuals); + +double ComputeResidualTime(std::vector vertexPos, double originTime, Hit lHitt); + +//* Gaussian kernel function +double GaussianKernel(double x, double bandwidth); + +//* Signal noise ratio +VtxCandidate CreateCandidate(const std::vector& vertex, double lowerLimit, double upperLimit, bool computeSNR); + +void VectorVertexPMT(std::vector vertex, int iPMT, double* dAngles ); + +std::vector ProjectPointToCylinder(std::vector point, double R, double H); \ No newline at end of file diff --git a/leaf/Likelihoods.cc b/leaf/Likelihoods.cc new file mode 100644 index 0000000..5cc79a2 --- /dev/null +++ b/leaf/Likelihoods.cc @@ -0,0 +1,743 @@ +#include "Likelihoods.hh" +#include "LeafInputs.hh" + +//Function to optimize, for MIGRAD +void MinuitLikelihood(int & /*nDim*/, double * /*gout*/, double &NLL, double par[], int /*flg*/) +{ + std::vector vertexPosition(4, 0.); // In centimeters + for (int i = 0; i < 4; i++) vertexPosition[i] = par[i]; + std::vector eventDirection(3, 0.); + for (int i = 0; i < 3; i++) eventDirection[i] = par[i + 10]; + // double dirThreshold = par[13]; + int nhits = par[5]; + double lowerLimit = par[6]; + double upperLimit = par[7]; + double directionality = par[9]; + // std::vector vertexDirection(3,0.); + double timeNLL = Likelihoods::Vertex_Time_NLL(fHitCollection, vertexPosition, nhits, lowerLimit, upperLimit, true, false, directionality); + // double angleNLL = Likelihoods::AngleNLL(fHitCollection, vertexPosition, fTrueDir); + NLL = timeNLL; +} + +//Function to optimize, for MIGRAD +void MinuitDirNLL(int& /*nDim*/, double* /*gout*/, double& DNLL, double par[], int /*flg*/) +{ + // Extract theta and phi from parameters + double theta = par[0]; + double phi = par[1]; + int nhits = static_cast(par[2]); + + // Ensure the vertexPosition vector is properly initialized + std::vector vertexPosition(4, 0.0); + for (int i = 0; i < 4; i++) vertexPosition[i] = par[i + 3]; + + // Ensure theta and phi are within valid ranges + if (theta < 0 || theta > TMath::Pi()) + { + DNLL = 1e10; // Assign a large NLL to invalid directions + return; + } + if (phi < -TMath::Pi() || phi > TMath::Pi()) + { + DNLL = 1e10; + return; + } + + // Calculate NLL using theta and phi + DNLL = Likelihoods::Dir_NLL(fHitCollection, vertexPosition, theta, phi, nhits); +} + +void MinuitJointNLL(int& nDim, double * gout, double & NLL, double par[], int flg) +{ + std::vector vertexPosition(4, 0.); + for (int i = 0; i < 4; i++) vertexPosition[i] = par[i]; + double theta = par[4]; + double phi = par[5]; + int nhits = par[7]; + double lowerLimit = par[8]; + double upperLimit = par[9]; + // double directionality = par[9]; + + // Ensure theta and phi are within valid ranges + if (theta < 0 || theta > TMath::Pi()) + { + NLL = 1e10; // Assign a large NLL to invalid directions + return; + } + if (phi < -TMath::Pi() || phi > TMath::Pi()) + { + NLL = 1e10; + return; + } + + double VtxNLL = Likelihoods::Vertex_Time_NLL(fHitCollection, vertexPosition, nhits, lowerLimit, upperLimit, true, false, 0); + double DirNLL = Likelihoods::Dir_NLL(fHitCollection, vertexPosition, theta, phi, nhits); + + NLL = VtxNLL * DirNLL; +} + +double Likelihoods::Vertex_Time_NLL(const HitCollection* lHitCol, std::vector vertexPosition, int nhits, double lowerLimit, double upperLimit, bool killEdges, bool scaleDR, int directionality) +{ + double NLL = 0; + + for (int ihit = 0; ihit < nhits; ihit++) + { + Hit lHit = lHitCol->At(ihit); + + int iPMT = lHit.PMT; + // double hitTime = (fTimeCorrection + lHit.T) / TimeDelta::ns; + PMTInfo lPMTInfo = (*fPMTList)[iPMT]; + + int pmtType = Astro_GetPMTType(iPMT); + // double distance = Astro_GetDistance(lPMTInfo.Position, vertexPosition); + + // double tof = distance / fLightSpeed; + // double residual = hitTime - tof - vertexPosition[3]; + double residual = ComputeResidualTime(vertexPosition, vertexPosition[3], lHitCol->At(ihit)); + + residual = ComputeResidualTime(vertexPosition, vertexPosition[3], lHit); + + // std::vector toHit = VectorToHitNorm(vertexPosition, lHit); + // if(dot(fTrueDir, toHit) < cos(90 * TMath::DegToRad())) + // if(dot(fTrueDir, toHit) < 0) + // { + // continue; //! Skip this hit if it's not in the same hemisphere + // } + + double proba = 0; + + bool condition = residual > lowerLimit && residual < upperLimit; + + if (condition) + { + proba = fSplineTimePDFQueue[pmtType]->Eval(residual); + +#ifdef VERBOSE_NLL + if (VERBOSE >= 3) + { + std::cout << "hit#" << ihit << ", hit time =" << hitTime << ", vertex time = " << vertexPosition[3] << std::endl; + std::cout << "Residual=" << residual << ", proba=" << proba << ", pmt type=" << pmtType << std::endl; + } +#endif + + if (scaleDR) + { + std::cout << "scaleDR with Lower Limit: " << lowerLimit << " Upper Limit: " << upperLimit << std::endl; + EventInfo fEventInfo = MakeEventInfo(lowerLimit, upperLimit, pmtType); + + // First, substract the DR from the PDF to keep only the signal: + double DR = fSplineTimePDFDarkRate[pmtType]->Eval(residual); + proba -= DR; + // Then, scale signal so that signal integral / DR integral = signalOverNoise + // To do so, we should scale signal so that integral = DR integral * signalOverNoise + // And of course, to rescale, we should divide by the signal integral + double Factor = fEventInfo.NoiseIntegral * fEventInfo.SignaloverNoise / fEventInfo.SignalIntegral; + proba *= Factor; + // And add again the DR + proba += DR; +#ifdef VERBOSE_NLL + if (VERBOSE >= 3) std::cout << "proba after scaling=" << proba << std::endl; +#endif + } + } + else if (killEdges) + { + if (residual >= upperLimit) + { + proba = fSplineTimePDFQueue[pmtType]->Eval(upperLimit); +#ifdef VERBOSE_NLL + if (VERBOSE >= 3 && pmtType == 1) + std::cout << "PMT type = " << pmtType << ", Upper limit = " << upperLimit << ", residual = " << residual << ", proba = " << proba << std::endl; +#endif + } + else if (residual <= lowerLimit) + { + // continue; + proba = fSplineTimePDFQueue[pmtType]->Eval(lowerLimit); +#ifdef VERBOSE_NLL + if (VERBOSE >= 3 && pmtType == 1) + std::cout << "PMT type = " << pmtType << ", Lower limit = " << lowerLimit << ", residual = " << residual << ", proba = " << proba << std::endl; +#endif + } + } + else proba = 0; + + if (proba < 0) + { + std::cout << "Error in PDF" << std::endl; + if (proba > -1e-1) proba = 0; // Since spline sometimes slightly goes below 0 due to interpolation + else + { + std::cout << "Error in " << residual << "ns where proba = " << proba << std::endl; + return 0; + } + } + if (proba == 0) proba = 1e-20; + NLL += -TMath::Log(proba); + } + + //* Directionality, not used RN + +// if (directionality != 0) +// { +// double NLLdir = 0; + +// NLLdir = this->FindNLLDirectionality(vertexPosition, nhits, VERBOSE, lowerLimit, upperLimit); + +// #ifdef VERBOSE_NLL +// if (VERBOSE >= 2) +// std::cout << "NLL = " << NLL << ", dir (L) = " << TMath::Exp(-NLLdir) << ", dir 2 = " << NLLdir << std::endl; +// #endif +// if (directionality == 1) +// NLL += NLLdir; +// else if (directionality == 2) +// NLL = NLLdir; +// } + + // timer.Stop(); + // std::cout << "NLL = " << NLL << std::endl; + return NLL; +} + +double Likelihoods::Vertex_Score(const HitCollection* lHitCol, std::vector vertexPosition, int nhits, double /*lowerLimit*/, double /*upperLimit*/, bool /*killEdges*/, bool /*scaleDR*/, int /*directionality*/) +{ + double NLL = 0; + + // std::cout << " Find NLL " << nhits << std::endl; + // std::cout << " First Hit " << fHitCollection->At(0).PMT << " " << vertexPosition.size() << std::endl; + for (int ihit = 0; ihit < nhits; ihit++) + { + // Hit lHit = fHitInfo[ihit]; + Hit lHit = lHitCol->At(ihit); + + // std::cout << " NLL Hit " << ihit << " " << lHit.PMT << std::endl; + int iPMT = lHit.PMT; + + // double hitTime = (fTimeCorrection + lHit.T) / TimeDelta::ns; + PMTInfo lPMTInfo = (*fPMTList)[iPMT]; + + int pmtType = Astro_GetPMTType(iPMT); + // double distance = Astro_GetDistance(lPMTInfo.Position, vertexPosition); + + // double tof = distance / fLightSpeed; + // double residual = hitTime - tof - vertexPosition[3]; + double residual = ComputeResidualTime(vertexPosition, vertexPosition[3], lHitCol->At(ihit)); +#ifdef KILLHALF + if (pmtType == 1) + { + double t = fRand->Uniform(0, 1); + if (t > 0.5) + continue; + } +#endif + + bool bCondition = (residual > fHitTimeLimitsNegative && residual < fHitTimeLimitsPositive) || (pmtType == 1 && !fLimit_mPMT); + + if (bCondition) + { + NLL++; +#ifdef VERBOSE_NLL + if (VERBOSE >= 3) + { + std::cout << "Residual=" << residual << ", pmt type=" << pmtType << std::endl; + std::cout << "hit#" << ihit << ", hit time =" << hitTime << ", vertex time = " << vertexPosition[3] << std::endl; + std::cout << "residual=" << residual << std::endl; + } +#endif + } + } + + return NLL != 0 ? -TMath::Log(NLL) : 1e15; +} + +double Likelihoods::FindNLL(const HitCollection* lHitCol, std::vector vertexPosition, int nhits, bool likelihood, int verbose, double lowerLimit, double upperLimit, bool killEdges, bool scaleDR, int directionality) +{ + double NLL = 0; + + for (int ihit = 0; ihit < nhits; ihit++) + { + Hit lHit = lHitCol->At(ihit); + + int iPMT = lHit.PMT; + // double hitTime = (fTimeCorrection + lHit.T) / TimeDelta::ns; + PMTInfo lPMTInfo = (*fPMTList)[iPMT]; + + int pmtType = Astro_GetPMTType(iPMT); + // double distance = Astro_GetDistance(lPMTInfo.Position, vertexPosition); + + // double tof = distance / fLightSpeed; + // double residual = hitTime - tof - vertexPosition[3]; + double residual = ComputeResidualTime(vertexPosition, vertexPosition[3], lHitCol->At(ihit)); + double proba; +#ifdef KILLHALF + if (pmtType == 1) + { + double t = fRand->Uniform(0, 1); + if (t > 0.5) + continue; + } +#endif + + if (likelihood) + { + bool condition; + + condition = residual > lowerLimit && residual < upperLimit; + if (condition) + { + proba = fSplineTimePDFQueue[pmtType]->Eval(residual); +#ifdef VERBOSE_NLL + if (VERBOSE >= 3) + { + std::cout << "hit#" << ihit << ", hit time =" << hitTime << ", vertex time = " << vertexPosition[3] << std::endl; + std::cout << "Residual=" << residual << ", proba=" << proba << ", pmt type=" << pmtType << std::endl; + } +#endif + if (scaleDR) + { + std::cout << "scaleDR with Lower Limit: " << lowerLimit << " Upper Limit: " << upperLimit << std::endl; + EventInfo fEventInfo = MakeEventInfo(lowerLimit, upperLimit, pmtType); + + // First, substract the DR from the PDF to keep only the signal: + double DR = fSplineTimePDFDarkRate[pmtType]->Eval(residual); + proba -= DR; + // Then, scale signal so that signal integral / DR integral = signalOverNoise + // To do so, we should scale signal so that integral = DR integral * signalOverNoise + // And of course, to rescale, we should divide by the signal integral + double Factor = fEventInfo.NoiseIntegral * fEventInfo.SignaloverNoise / fEventInfo.SignalIntegral; + proba *= Factor; + // And add again the DR + proba += DR; + if (VERBOSE >= 3) + std::cout << "proba after scaling=" << proba << std::endl; + } + } + else if (killEdges && residual >= upperLimit) + { + proba = fSplineTimePDFQueue[pmtType]->Eval(upperLimit); +#ifdef VERBOSE_NLL + if (VERBOSE >= 3 && pmtType == 1) + std::cout << "PMT type = " << pmtType << ", Upper limit = " << upperLimit << ", residual = " << residual << ", proba = " << proba << std::endl; +#endif + } + else if (killEdges && residual <= lowerLimit) + { + proba = fSplineTimePDFQueue[pmtType]->Eval(lowerLimit); +#ifdef VERBOSE_NLL + if (VERBOSE >= 3 && pmtType == 1) + std::cout << "PMT type = " << pmtType << ", Lower limit = " << lowerLimit << ", residual = " << residual << ", proba = " << proba << std::endl; +#endif + } + else + { + proba = 0; + } + } + else + { + bool condition; + if (pmtType == 0 || (pmtType == 1 && fLimit_mPMT)) condition = residual > fHitTimeLimitsNegative && residual < fHitTimeLimitsPositive; + else condition = true; + if (condition) + { + NLL++; //= fSplineTimePDFQueue[pmtType]->Eval(residual); +#ifdef VERBOSE_NLL + if (VERBOSE >= 3) + { + std::cout << "Residual=" << residual << ", pmt type=" << pmtType << std::endl; + std::cout << "hit#" << ihit << ", hit time =" << hitTime << ", vertex time = " << vertexPosition[3] << ", distance PMT vs vertex = " << distance << ", tof=" << tof << std::endl; + std::cout << "residual=" << residual << ", proba=" << proba << std::endl; + } +#endif + } + else if (killEdges) continue; + else proba = 0; + } + + if (likelihood) + { + if (proba < 0) + { + std::cout << "Error in PDF" << std::endl; + if (proba > -1e-1) proba = 0; // Since spline sometimes slightly goes below 0 due to interpolation + else + { + std::cout << "Error in " << residual << "ns where proba = " << proba << std::endl; + return 0; + } + } + if (proba == 0) proba = 1e-20; + NLL += -TMath::Log(proba); + } + } + // std::cout << "FindNLL: Time it took for the loop = " << timer.RealTime() << std::endl; + // cout<<"NLL="<=2) cout<<"NLL="<38) + { + std::cout << "Got hit! " << theta << std::endl; + DNLL += -1; + } + + + if (VERBOSE >= 2) std::cout << "Direction NLL = " << DNLL << std::endl; + } + + return DNLL; + +} + +//Uses the PDF to find the DirNLL, used mainly for the MIGRAD optimization +double Likelihoods::Dir_NLL(const HitCollection* lHitCol, const std::vector& vertexPosition, double theta_track, double phi_track, int nhits) { + double DNLL = 0; + + // Ensure theta and phi are within valid ranges + // theta in [0, pi], phi in [-pi, pi] + if (theta_track < 0 || theta_track > TMath::Pi()) + { + std::cerr << "Error: theta_track is out of range [0, π]." << std::endl; + return std::numeric_limits::infinity(); + } + if (phi_track < -TMath::Pi() || phi_track > TMath::Pi()) + { + std::cerr << "Error: phi_track is out of range [-π, π]." << std::endl; + return std::numeric_limits::infinity(); + } + + // Precompute sine and cosine of theta_track and phi_track + // double sin_theta_track = sin(theta_track); + // double cos_theta_track = cos(theta_track); + // double sin_phi_track = sin(phi_track); + // double cos_phi_track = cos(phi_track); + /* 1 + double a = 0.556805; + double b = 1.06697; + double c = -0.902084; + double d = -0.0707118;*/ + + // 2 + // double a = -0.000440808; + // double b = 1.71313; + // double c = -2.01675; + // double d = 1.81755; + + // Direction vector components from theta and phi + // double vertexDirNormalized[3]; + // vertexDirNormalized[0] = sin_theta_track * cos_phi_track; + // vertexDirNormalized[1] = sin_theta_track * sin_phi_track; + // vertexDirNormalized[2] = cos_theta_track; + + for (int ihit = 0; ihit < nhits; ihit++) { + int trueNHits = lHitCol->Size(); + if(ihit >= trueNHits) + { + std::cout << "Error: ihit is out of range in direction NLL computation" << std::endl; + return std::numeric_limits::infinity(); + } + if(vertexPosition.size() != 4) + { + std::cout << "Error: vertexPosition is not of size 4 in direction NLL computation" << std::endl; + return std::numeric_limits::infinity(); + } + + Hit lHit = lHitCol->At(ihit); + int iPMT = lHit.PMT; + + // double hitTime = (fTimeCorrection + lHit.T) / TimeDelta::ns; + PMTInfo lPMTInfo = (*fPMTList)[iPMT]; + + int pmtType = Astro_GetPMTType(iPMT); + // double distance = Astro_GetDistance(lPMTInfo.Position,vertexPosition); + + // double tof = distance / fLightSpeed; + // double residual = hitTime - tof - vertexPosition[3]; + + double residual = ComputeResidualTime(vertexPosition, vertexPosition[3], lHitCol->At(ihit)); + + bool DirCondition = (residual >= fDirectionPDF_minResidual && residual <= fDirectionPDF_maxResidual) || (pmtType == 1 && !fLimit_mPMT); + + + if(DirCondition || DirTakeAll) + { + int iPMT = lHit.PMT; + PMTInfo lPMTInfo = (*fPMTList)[iPMT]; + + // Ensure lPMTInfo.Position and lPMTInfo.Orientation have the correct size + if (sizeof(lPMTInfo.Position) / sizeof(lPMTInfo.Position[0]) < 3 || + sizeof(lPMTInfo.Orientation) / sizeof(lPMTInfo.Orientation[0]) < 3) { + std::cerr << "Error: lPMTInfo.Position or lPMTInfo.Orientation is not of size 3." << std::endl; + continue; + } + + double dirVector[3]; + dirVector[0] = lPMTInfo.Position[0] - vertexPosition[0]; + dirVector[1] = lPMTInfo.Position[1] - vertexPosition[1]; + dirVector[2] = lPMTInfo.Position[2] - vertexPosition[2]; + + double dirNorm = sqrt(dirVector[0]*dirVector[0] + dirVector[1]*dirVector[1] + dirVector[2]*dirVector[2]); + if (dirNorm == 0) + { + std::cerr << "Error: dirNorm is zero, cannot normalize direction vector." << std::endl; + continue; + } + dirVector[0] /= dirNorm; + dirVector[1] /= dirNorm; + dirVector[2] /= dirNorm; + + // Convert dirVector to spherical coordinates (theta_pmt, phi_pmt)fDirectionPDF + double theta_pmt = acos(dirVector[2]); // theta in [0, π] + double phi_pmt = atan2(dirVector[1], dirVector[0]); // phi in [-π, π] + + // Calculate the angular difference between the two directions + // double delta_theta = theta_pmt - theta_track; //* could be useful + double delta_phi = phi_pmt - phi_track; + + // Ensure delta_phi is within [-π, π] + if (delta_phi > TMath::Pi()) delta_phi -= 2 * TMath::Pi(); + else if (delta_phi < -TMath::Pi()) delta_phi += 2 * TMath::Pi(); + + // Calculate the angle between the two directions using the spherical law of cosines + double cos_relative_angle = sin(theta_track) * sin(theta_pmt) * cos(delta_phi) + + cos(theta_track) * cos(theta_pmt); + + //cos_relative_angle = std::max(-1.0, std::min(1.0, cos_relative_angle)); + + double relative_angle = acos(cos_relative_angle); // Angle in radians + + // Convert angle to degrees if your PDF is defined in degrees + double theta = relative_angle * 180.0 / TMath::Pi(); + + // Evaluate your PDF at theta + // std::cout << "before pdf : " << ihit << std::endl; + // double proba = fDirectionPDF->Eval(cos_relative_angle)/*sin(theta)*/; //* if cos spline + // double proba = fDirectionPDF->Eval(cos_relative_angle); //* Taha + double proba = fHitAnglePDF->Eval(theta); //* Nicolas + // std::cout << "after pdf : " << ihit << std::endl; + + // Correct for PMT orientation if necessary + double PMTOrientation[3]; + PMTOrientation[0] = lPMTInfo.Orientation[0]; + PMTOrientation[1] = lPMTInfo.Orientation[1]; + PMTOrientation[2] = lPMTInfo.Orientation[2]; + + // Normalize PMT orientation vector + double pmtOrientNorm = sqrt(PMTOrientation[0]*PMTOrientation[0] + + PMTOrientation[1]*PMTOrientation[1] + + PMTOrientation[2]*PMTOrientation[2]); + + PMTOrientation[0] /= pmtOrientNorm; + PMTOrientation[1] /= pmtOrientNorm; + PMTOrientation[2] /= pmtOrientNorm; + + double cosPMThitAngle = dirVector[0]*PMTOrientation[0] + + dirVector[1]*PMTOrientation[1] + + dirVector[2]*PMTOrientation[2]; + // double PMThitAngle = acos(-cosPMThitAngle); + cosPMThitAngle *= -1; + // double cos2PMThitAngle = fabs(-2 * cosPMThitAngle * cosPMThitAngle + 1); + // double Ftheta = a + b*cosPMThitAngle + c*(pow(cosPMThitAngle,2)) + d*(pow(cosPMThitAngle,3)); + + //proba *= hitCharge; + + // Avoid log(0) by ensuring proba is positive + if (proba <= 0) + { + proba = 1e-20; + if (VERBOSE >= 1) std::cout << "Warning: Probability is zero or negative at hit #" << ihit << std::endl; + } + //std::cout << "Proba Corrected by " << PMThitAngle * 180 / M_PI << " degrees with a sincos " << sin(PMThitAngle) << ", " << cosPMThitAngle << std::endl; + // Accumulate Negative Log-Likelihood + DNLL += -TMath::Log(proba);//*(cosPMThitAngle/Ftheta)*(sin(PMThitAngle)); + + if (VERBOSE >= 3) std::cout << "Hit #" << ihit << ", Theta = " << theta << " degrees, Probability = " << proba << std::endl; + } + + if (VERBOSE >= 2) std::cout << "Total Direction NLL = " << DNLL << std::endl; + } + return DNLL; +} + +double Likelihoods::AngleNLL(const HitCollection* lHitCol, std::vector vertexPosition, std::vector vertexDirection) +{ + double NLL = 0; + + int nhits = lHitCol->Size(); + + for (int ihit = 0; ihit < nhits; ihit++) + { + //* compute angle + Hit lHit = lHitCol->At(ihit); + int iPMT = lHit.PMT; + PMTInfo lPMTInfo = (*fPMTList)[iPMT]; + std::vector toPMT = std::vector(3); + for(int j = 0; j < 3; j++) toPMT[j] = lPMTInfo.Position[j] - vertexPosition[j]; + Normalize(toPMT); + Normalize(vertexDirection); + double dotP = dot(toPMT, vertexDirection); + dotP = std::max(-1.0, std::min(1.0, dotP)); // Ensure dot product is within valid range for acos + double angle = std::acos(dotP) * 180 / TMath::Pi(); + // std::cout<<"angle : " << angle << " on pdf : " << fHitAnglePDF->Eval(angle) << std::endl; + // NLL += -TMath::Log(std::max(1e-20, fDirectionPDF->Eval(dotP))); //? Taha + NLL += -TMath::Log(std::max(1e-20, fHitAnglePDF->Eval(angle))); //? Nicolas + } + + return NLL; +} + + +double Likelihoods::FindNLLDirectionality(const HitCollection* lHitCol, std::vector vVtxPos, int nhits, int verbose, double /*lowerLimit*/, double /*upperLimit*/) +{ + + double NLL = 0; + // std::vector vDirection(2, 0.); // Return phi and theta. + + // // Interpolate isn't compatible with multi-thread, need to use mutex which lead to long deadtime. + // // Copy TGraph2D + // mtx.lock(); + // static thread_local TGraph2D gPMTDirectionality_2D_local_0 = TGraph2D(*gPMTDirectionality_2D[MiniPMT][0]); + // static thread_local TGraph2D gPMTDirectionality_2D_local_1 = TGraph2D(*gPMTDirectionality_2D[MiniPMT][1]); + // static thread_local TGraph2D gPMTDirectionality_2D_local_2 = TGraph2D(*gPMTDirectionality_2D[MiniPMT][2]); + // mtx.unlock(); + + // TGraph2D *tDirectionality[3]; + // tDirectionality[0] = &gPMTDirectionality_2D_local_0; + // tDirectionality[1] = &gPMTDirectionality_2D_local_1; + // tDirectionality[2] = &gPMTDirectionality_2D_local_2; + + // for (int ihit = 0; ihit < nhits; ihit++) + // { + // // Hit lHit = fHitInfo[ihit]; + // Hit lHit = LeafInputs::fHitCollection->At(ihit); + + // int iPMT = lHit.PMT; + + // PMTInfo lPMTInfo = (*fPMTList)[iPMT]; + // int pmtType = Astro_GetPMTType(iPMT); + + // if (pmtType == 0) + // continue; + // bool condition = true; + + // if (condition) + // { + // double vPMTVtx[4]; + // this->VectorVertexPMT(vVtxPos, iPMT, vPMTVtx); + + // double dPhi = vPMTVtx[0]; + // double dTheta = vPMTVtx[1]; + // double dDist = vPMTVtx[2]; + + // int pmtGroup = lPMTInfo.mPMT_Group; + + // double proba = 0; + // proba = tDirectionality[pmtGroup]->Interpolate(dPhi, dTheta); + + // double dDistCorr = fDistResponsePMT[pmtType]->Eval(dDist); + + // proba *= dDistCorr; + + // if (proba == 0) + // { + // proba = 1e-20; + // if (verbose) + // std::cout << "We are at proba = 0, theta = " << dTheta << ", proba used = " << proba << std::endl; + // } + // NLL += -TMath::Log(proba); + + // if (VERBOSE >= 3) + // { + // // std::cout<<"PMT type="< + struct CompareByNLL { + bool operator()(T const &a, T const &b) const { + return a.NLL < b.NLL; + } + }; + + //* Calculate likelihood without using PDF, but just using hits within a given timing window (hits "in-time"). + static double Vertex_Score(const HitCollection* lHitCol, std::vector vertexPosition, int nhits, double lowerLimit, double upperLimit, bool killEdges, bool scaleDR, int directionality); + + //* Contain the two functions above in one function, where usage of PDF or not can be set through the flag: likelihood = true/false + static double FindNLL(const HitCollection* lHitCol, std::vector vertexPosition,int nhits, bool likelihood, int verbose, double lowerLimit, double upperLimit, bool killEdges=false, bool scaleDR=false, int directionality=false); + + static double Vertex_Time_NLL(const HitCollection* lHitCol, std::vector vertexPosition, int nhits, double lowerLimit, double upperLimit, bool killEdges, bool scaleDR, int directionality); + + static double Dir_NLL(const HitCollection* lHitCol, const std::vector& vertexPosition, const double theta_track, const double phi_track, int nhits); + + static double ComputeDirNLL_NoPDF(const HitCollection* lHitCol, const std::vector& vertexPosition, const std::vector& vertexDirection, int nhits); + + static double AngleNLL(const HitCollection* lHitCol, std::vector vertexPosition, std::vector vertexDirection); + + static double FindNLLDirectionality(const HitCollection* lHitCol, std::vector vertexPosition, int nhits, int verbose, double lowerLimit, double upperLimit); + + static double GoodnessOfFit(const HitCollection* lHitCol, std::vector vertexPosition, int nhits, double lowerLimit, double upperLimit, bool killEdges, bool scaleDR, int directionality); + +}; \ No newline at end of file