From b2c650b7980f7338622f27ca92203e0254ccaa51 Mon Sep 17 00:00:00 2001 From: flemmons <85371521+flemmons@users.noreply.github.com> Date: Fri, 16 Feb 2024 16:16:20 -0700 Subject: [PATCH] Update to Vertex Reconstruction (#226) * Update to Vertex Reconstruction, including new FineGrid search and mean time selection * False Commit to enable checkout; please ignore. * Edited recent vertex reconstruction updates for better mesh with TA code Removed redundant GenerateFineGrid function in VtxSeedGenerator Removed experimental penalty term calculation in FoMCalculator * Fixed small typo in VtxSeedGenerator.cpp * Fixes to VtxSeedFineGrid. Added new PDF charge fitter to FoMCalculator, and handling in VtxExtendedVertexFinder and MinuitOptimizer. This change is experimental, but should only trigger with appropriate config variables. * bug fixes to recent update * Fixes to VtxExtendedVertexFinder, and further progress in experimental charge fit. * Incremental update on pdf charge fit. Updated to chi2 definition. * Charge fit update for LikelihoodCheck tool, along with a debug ability to vary 2D angle instead of position * Further debug update of charge fit addition to FoMCalculator * Bringing up to date with new charge fitter fixes. Added VGCheck-pdfcharge.root, which is the file containing pdf plot for use within new charge fitter. --- DataModel/FoMCalculator.cpp | 1181 +++++++++-------- DataModel/FoMCalculator.h | 6 +- DataModel/MinuitOptimizer.cpp | 149 +++ DataModel/MinuitOptimizer.h | 1 + UserTools/Factory/Factory.cpp | 2 + .../LikelihoodFitterCheck.cpp | 112 +- .../LikelihoodFitterCheck.h | 10 + UserTools/LikelihoodFitterCheck/README.md | 0 UserTools/MeanTimeCheck/MeanTimeCheck.cpp | 185 +++ UserTools/MeanTimeCheck/MeanTimeCheck.h | 70 + UserTools/MeanTimeCheck/README.md | 20 + UserTools/Unity.h | 2 + UserTools/VertexGeometryCheck/README.md | 0 .../VertexGeometryCheck.cpp | 10 +- .../VertexGeometryCheck/VertexGeometryCheck.h | 1 + UserTools/VtxExtendedVertexFinder/README.md | 0 .../VtxExtendedVertexFinder.cpp | 32 +- .../VtxExtendedVertexFinder.h | 10 +- UserTools/VtxSeedFineGrid/README.md | 20 + UserTools/VtxSeedFineGrid/VtxSeedFineGrid.cpp | 412 ++++++ UserTools/VtxSeedFineGrid/VtxSeedFineGrid.h | 85 ++ UserTools/VtxSeedGenerator/README.md | 0 .../VtxSeedGenerator/VtxSeedGenerator.cpp | 8 +- UserTools/VtxSeedGenerator/VtxSeedGenerator.h | 10 + VGCheck-pdfcharge.root | Bin 0 -> 19920 bytes 25 files changed, 1767 insertions(+), 559 deletions(-) mode change 100644 => 100755 DataModel/FoMCalculator.cpp mode change 100644 => 100755 DataModel/FoMCalculator.h mode change 100644 => 100755 DataModel/MinuitOptimizer.cpp mode change 100644 => 100755 DataModel/MinuitOptimizer.h mode change 100644 => 100755 UserTools/LikelihoodFitterCheck/LikelihoodFitterCheck.cpp mode change 100644 => 100755 UserTools/LikelihoodFitterCheck/LikelihoodFitterCheck.h mode change 100644 => 100755 UserTools/LikelihoodFitterCheck/README.md create mode 100644 UserTools/MeanTimeCheck/MeanTimeCheck.cpp create mode 100644 UserTools/MeanTimeCheck/MeanTimeCheck.h create mode 100644 UserTools/MeanTimeCheck/README.md mode change 100644 => 100755 UserTools/VertexGeometryCheck/README.md mode change 100644 => 100755 UserTools/VertexGeometryCheck/VertexGeometryCheck.cpp mode change 100644 => 100755 UserTools/VertexGeometryCheck/VertexGeometryCheck.h mode change 100644 => 100755 UserTools/VtxExtendedVertexFinder/README.md mode change 100644 => 100755 UserTools/VtxExtendedVertexFinder/VtxExtendedVertexFinder.cpp mode change 100644 => 100755 UserTools/VtxExtendedVertexFinder/VtxExtendedVertexFinder.h create mode 100755 UserTools/VtxSeedFineGrid/README.md create mode 100755 UserTools/VtxSeedFineGrid/VtxSeedFineGrid.cpp create mode 100755 UserTools/VtxSeedFineGrid/VtxSeedFineGrid.h mode change 100644 => 100755 UserTools/VtxSeedGenerator/README.md mode change 100644 => 100755 UserTools/VtxSeedGenerator/VtxSeedGenerator.cpp mode change 100644 => 100755 UserTools/VtxSeedGenerator/VtxSeedGenerator.h create mode 100644 VGCheck-pdfcharge.root diff --git a/DataModel/FoMCalculator.cpp b/DataModel/FoMCalculator.cpp old mode 100644 new mode 100755 index d2f859a64..57d77b59f --- a/DataModel/FoMCalculator.cpp +++ b/DataModel/FoMCalculator.cpp @@ -1,529 +1,652 @@ -#include "FoMCalculator.h" - -//Constructor -FoMCalculator::FoMCalculator() { - fVtxGeo = 0; - fBaseFOM = 100.0; - - // fitting parameters ported from vertexFinder (FIXME) - fTimeFitWeight = 0.50; // nominal time weight - fConeFitWeight = 0.50; // nominal cone weight - - // default Mean time calculator type - fMeanTimeCalculatorType = 0; -} - -//Destructor -FoMCalculator::~FoMCalculator() { - fVtxGeo = 0; -} - - -void FoMCalculator::LoadVertexGeometry(VertexGeometry* vtxgeo) { - this->fVtxGeo = vtxgeo; -} - - -void FoMCalculator::TimePropertiesLnL(double vtxTime, double& vtxFOM) -{ - // internal variables - // ================== - double weight = 0.0; - double delta = 0.0; // time residual of each hit - double sigma = 0.0; // time resolution of each hit - double A = 0.0; // normalisation of first Gaussian - int type; // Digit type (LAPPD or PMT) - - double Preal = 0.0; // probability of real hit - double P = 0.0; // probability of hit - - double chi2 = 0.0; // log-likelihood: chi2 = -2.0*log(L) - double ndof = 0.0; // total number of hits - double fom = -9999.; // figure of merit - - // add noise to model - // ================== - double Pnoise; - - //FIXME: We need an implementation of noise models for PMTs and LAPPDs - //double nFilterDigits = this->fVtxGeo->GetNFilterDigits(); - //double fTimeFitNoiseRate = 0.02; // hits/ns [0.40 for electrons, 0.02 for muons] - //Pnoise = fTimeFitNoiseRate/nFilterDigits; - - // loop over digits - // ================ - for( int idigit=0; idigitfVtxGeo->GetNDigits(); idigit++ ){ - int detType = this->fVtxGeo->GetDigitType(idigit); - delta = this->fVtxGeo->GetDelta(idigit) - vtxTime; - sigma = this->fVtxGeo->GetDeltaSigma(idigit); - type = this->fVtxGeo->GetDigitType(idigit); - if (type == RecoDigit::PMT8inch){ //PMT8Inch - sigma = 1.5*sigma; - Pnoise = 1e-8; //FIXME; Need implementation of noise model - } - if (type == RecoDigit::lappd_v0){ //lappd - sigma = 1.2*sigma; - Pnoise = 1e-8; //FIXME; Need implementation of noise model - } - A = 1.0 / ( 2.0*sigma*sqrt(0.5*TMath::Pi()) ); //normalisation constant - Preal = A*exp(-(delta*delta)/(2.0*sigma*sigma)); - P = (1.0-Pnoise)*Preal + Pnoise; - chi2 += -2.0*log(P); - ndof += 1.0; - } - - // calculate figure of merit - if( ndof>0.0 ){ - fom = fBaseFOM - 5.0*chi2/ndof; - } - - // return figure of merit - // ====================== - vtxFOM = fom; - return; -} - - - - -void FoMCalculator::ConePropertiesFoM(double coneEdge, double& coneFOM) -{ - // calculate figure of merit - // ========================= - double coneEdgeLow = 21.0; // cone edge (low side) - double coneEdgeHigh = 3.0; // cone edge (high side) [muons: 3.0, electrons: 7.0] - double deltaAngle = 0.0; - double digitCharge = 0.0; - double coneCharge = 0.0; - double allCharge = 0.0; - - double fom = -9999.; - - for( int idigit=0; idigitfVtxGeo->GetNDigits(); idigit++ ){ - if( this->fVtxGeo->IsFiltered(idigit) && this->fVtxGeo->GetDigitType(idigit) == RecoDigit::PMT8inch){ - deltaAngle = this->fVtxGeo->GetAngle(idigit) - coneEdge; - digitCharge = this->fVtxGeo->GetDigitQ(idigit); - - if( deltaAngle<=0.0 ){ - coneCharge += digitCharge*( 0.75 + 0.25/( 1.0 + (deltaAngle*deltaAngle)/(coneEdgeLow*coneEdgeLow) ) ); - } - else{ - coneCharge += digitCharge*( 0.00 + 1.00/( 1.0 + (deltaAngle*deltaAngle)/(coneEdgeHigh*coneEdgeHigh) ) ); - } - - allCharge += digitCharge; - } - } - - if( allCharge>0.0 ){ - fom = fBaseFOM*coneCharge/allCharge; - } - - // return figure of merit - // ====================== - coneFOM = fom; - return; -} - - - - -//Given the position of the point vertex (x, y, z) and n digits, calculate the mean expected vertex time -double FoMCalculator::FindSimpleTimeProperties(double myConeEdge) { - double meanTime = 0.0; - // weighted average - if(fMeanTimeCalculatorType == 0) { - // calculate mean and rms of hits inside cone - // ========================================== - double Swx = 0.0; - double Sw = 0.0; - - double delta = 0.0; - double sigma = 0.0; - double weight = 0.0; - double deweight = 0.0; - double deltaAngle = 0.0; - - double myConeEdgeSigma = 7.0; // [degrees] - - for( int idigit=0; idigitfVtxGeo->GetNDigits(); idigit++ ){ - int detType = this->fVtxGeo->GetDigitType(idigit); - if( this->fVtxGeo->IsFiltered(idigit) ){ - delta = this->fVtxGeo->GetDelta(idigit); - sigma = this->fVtxGeo->GetDeltaSigma(idigit); - weight = 1.0/(sigma*sigma); - // profile in angle - deltaAngle = this->fVtxGeo->GetAngle(idigit) - myConeEdge; - // deweight hits outside cone - if( deltaAngle<=0.0 ){ - deweight = 1.0; - } - else{ - deweight = 1.0/(1.0+(deltaAngle*deltaAngle)/(myConeEdgeSigma*myConeEdgeSigma)); - } - Swx += deweight*weight*delta; //delta is expected vertex time - Sw += deweight*weight; - } - } - if( Sw>0.0 ){ - meanTime = Swx*1.0/Sw; - } - } - - // most probable time - else if(fMeanTimeCalculatorType == 1) { - double sigma = 0.0; - double deltaAngle = 0.0; - double weight = 0.0; - double deweight = 0.0; - double myConeEdgeSigma = 7.0; // [degrees] - vector deltaTime1; - vector deltaTime2; - vector TimeWeight; - - for( int idigit=0; idigitGetNDigits(); idigit++ ){ - if(fVtxGeo->IsFiltered(idigit)){ - deltaTime1.push_back(fVtxGeo->GetDelta(idigit)); - deltaTime2.push_back(fVtxGeo->GetDelta(idigit)); - sigma = fVtxGeo->GetDeltaSigma(idigit); - weight = 1.0/(sigma*sigma); - deltaAngle = fVtxGeo->GetAngle(idigit) - myConeEdge; - if( deltaAngle<=0.0 ){ - deweight = 1.0; - } - else{ - deweight = 1.0/(1.0+(deltaAngle*deltaAngle)/(myConeEdgeSigma*myConeEdgeSigma)); - } - TimeWeight.push_back(deweight*weight); - } - } - int n = deltaTime1.size(); - std::sort(deltaTime1.begin(),deltaTime1.end()); - double timeMin = deltaTime1.at(int((n-1)*0.05)); // 5% of the total entries - double timeMax = deltaTime1.at(int((n-1)*0.90)); // 90% of the total entries - int nbins = int(n/5); - TH1D *hDeltaTime = new TH1D("hDeltaTime", "hDeltaTime", nbins, timeMin, timeMax); - for(int i=0; iFill(deltaTime2.at(i), TimeWeight.at(i)); - hDeltaTime->Fill(deltaTime2.at(i)); - } - meanTime = hDeltaTime->GetBinCenter(hDeltaTime->GetMaximumBin()); - delete hDeltaTime; hDeltaTime = 0; - } - - else std::cout<<"FoMCalculator Error: Wrong type of Mean time calculator! "<fVtxGeo->CalcPointResiduals(vtxX, vtxY, vtxZ, 0.0, 0.0, 0.0, 0.0); //calculate expected point vertex time for each digit - - // calculate figure of merit - // ========================= - this->TimePropertiesLnL(vtxTime, vtxFOM); - - // calculate overall figure of merit - // ================================= - fom = vtxFOM; - // truncate - if( fom<-9999. ) fom = -9999.; - - return; -} - - - - - -void FoMCalculator::PointDirectionChi2(double vtxX, double vtxY, - double vtxZ, double dirX, double dirY, double dirZ, - double coneAngle, double& fom) -{ - // figure of merit - // =============== - double coneFOM = -9999.; - - // calculate residuals - // =================== - this->fVtxGeo->CalcPointResiduals(vtxX, vtxY, vtxZ, 0.0, - dirX, dirY, dirZ); //load expected vertex time for each digit - - // calculate figure of merit - // ========================= - this->ConePropertiesFoM(coneAngle, coneFOM); - - // calculate overall figure of merit - // ================================= - fom = coneFOM; - - // truncate - if( fom<-9999. ) fom = -9999.; - - return; -} - -void FoMCalculator::PointVertexChi2(double vtxX, double vtxY, double vtxZ, - double dirX, double dirY, double dirZ, - double coneAngle, double vtxTime, double& fom) -{ - // figure of merit - // =============== - double vtxFOM = -9999.; - - // calculate residuals - // =================== - this->fVtxGeo->CalcPointResiduals(vtxX, vtxY, vtxZ, 0.0, - dirX, dirY, dirZ); //calculate expected vertex time for each digit - // calculate figure of merit - // ========================= - double timeFOM = -9999.; - double coneFOM = -9999.; - - this->ConePropertiesFoM(coneAngle,coneFOM); - this->TimePropertiesLnL(vtxTime, timeFOM); - - double fTimeFitWeight = this->fTimeFitWeight; - double fConeFitWeight = this->fConeFitWeight; - vtxFOM = (fTimeFitWeight*timeFOM+fConeFitWeight*coneFOM)/(fTimeFitWeight+fConeFitWeight); - - // calculate overall figure of merit - // ================================= - fom = vtxFOM; - // truncate - if( fom<-9999. ) fom = -9999.; - - return; -} - -void FoMCalculator::ExtendedVertexChi2(double vtxX, double vtxY, double vtxZ, double dirX, double dirY, double dirZ, double coneAngle, double vtxTime, double& fom) -{ - // figure of merit - // =============== - double vtxFOM = -9999.; - double timeFOM = -9999.; - double coneFOM = -9999.; - - // calculate residuals - // =================== - this->fVtxGeo->CalcExtendedResiduals(vtxX,vtxY,vtxZ,0.0,dirX,dirY,dirZ); - - // calculate figure of merit - // ========================= - - this->ConePropertiesFoM(coneAngle,coneFOM); - this->TimePropertiesLnL(vtxTime, timeFOM); - - double fTimeFitWeight = this->fTimeFitWeight; - double fConeFitWeight = this->fConeFitWeight; - vtxFOM = (fTimeFitWeight*timeFOM+fConeFitWeight*coneFOM)/(fTimeFitWeight+fConeFitWeight); - - // calculate overall figure of merit - // ================================= - fom = vtxFOM; - - // truncate - if( fom<-9999. ) fom = -9999.; - - return; -} - -//KEPT FOR HISTORY, BUT FITTER IS CURRENTLY NOT WORKING -//void FoMCalculator::CorrectedVertexChi2(double vtxX, double vtxY, double vtxZ, double dirX, double dirY, double dirZ, double& vtxAngle, double& vtxTime, double& fom) -//{ -// // figure of merit -// // =============== -// double vtxFOM = 0.0; -// double timeFOM = 0.0; -// double coneFOM = 0.0; -// double penaltyFOM = 0.0; -// double fixPositionFOM = 0.0; -// double fixDirectionFOM = 0.0; -// -// // calculate residuals -// // =================== -// this->fVtxGeo->CalcExtendedResiduals(vtxX,vtxY,vtxZ,0.0,dirX,dirY,dirZ); -// -// // calculate figure of merit -// // ========================= -// -// this->ConePropertiesFoM(vtxAngle,coneFOM); -// fom = coneFOM; -// -// // truncate -// if( fom<-999.999*fBaseFOM ) fom = -999.999*fBaseFOM; -// -// return; -//} -// -//void FoMCalculator::ConePropertiesLnL(double coneParam0, double coneParam1, double coneParam2, double& coneAngle, double& coneFOM) -//{ -// -// // nuisance parameters -// // =================== -// double alpha = coneParam0; // track length parameter = 0.25 -// double alpha0 = coneParam1; // track length parameter = 0.5 -// double beta = coneParam2; // particle ID: 0.0[electron]->1.0[muon] = 0.75 -// -// // internal variables -// // ================== -// double deltaAngle = 0.0; // -// double sigmaAngle = 7.0; //Cherenkov Angle resolution -// double Angle0 = Parameters::CherenkovAngle(); //Cherenkov Angle: 42 degree -// double deltaAngle0 = Angle0*alpha0; //? -// -// double digitQ = 0.0; -// double sigmaQmin = 1.0; -// double sigmaQmax = 10.0; -// double sigmaQ = 0.0; -// -// double A = 0.0; -// -// double PconeA = 0.0; -// double PconeB = 0.0; -// double Pmu = 0.0; -// double Pel = 0.0; -// -// double Pcharge = 0.0; -// double Pangle = 0.0; -// double P = 0.0; -// -// double chi2 = 0.0; -// double ndof = 0.0; -// -// double angle = 46.0; -// double fom = 0.0; -// -// // hard-coded parameters: 200 kton (100 kton) -// // ========================================== -// double lambdaMuShort = 0.5; // 0.5; -// double lambdaMuLong = 5.0; // 15.0; -// double alphaMu = 1.0; // 4.5; -// -// double lambdaElShort = 1.0; // 2.5; -// double lambdaElLong = 7.5; // 15.0; -// double alphaEl = 6.0; // 3.5; -// -// // numerical integrals -// // =================== -// fSconeA = 21.9938; -// fSconeB = 0.0000; -// -// // Number of P.E. angular distribution -// // inside cone -// int nbinsInside = 420; //divide the angle range by 420 bins -// for( int n=0; nfVtxGeo->GetNDigits(); idigit++ ){ -// -// if( this->fVtxGeo->IsFiltered(idigit) ){ -// digitQ = this->fVtxGeo->GetDigitQ(idigit); -// deltaAngle = this->fVtxGeo->GetAngle(idigit) - Angle0; -// -// // pulse height distribution -// // ========================= -// if( deltaAngle<=0 ){ //inside Cone -// sigmaQ = sigmaQmax; -// } -// else{ //outside Cone -// sigmaQ = sigmaQmin + (sigmaQmax-sigmaQmin)/(1.0+(deltaAngle*deltaAngle)/(sigmaAngle*sigmaAngle)); -// } -// -// A = 1.0/(log(2.0)+0.5*TMath::Pi()*sigmaQ); -// -// if( digitQ<1.0 ){ -// Pcharge = 2.0*A*digitQ/(1.0+digitQ*digitQ); -// } -// else{ -// Pcharge = A/(1.0+(digitQ-1.0)*(digitQ-1.0)/(sigmaQ*sigmaQ)); -// } -// -// // angular distribution -// // ==================== -// A = 1.0/( alpha*fSconeA + (1.0-alpha)*fSconeB -// + beta*fSmu + (1.0-beta)*fSel ); // numerical integrals -// -// if( deltaAngle<=0 ){ -// -// // pdfs inside cone: -// PconeA = 1.4944765*sin( (Angle0+deltaAngle)*(TMath::Pi()/180.0) ); -// PconeB = 1.4944765*sin( (Angle0+deltaAngle)*(TMath::Pi()/180.0) ) -// *( 1.0/(1.0+(deltaAngle*deltaAngle)/(deltaAngle0*deltaAngle0)) ); -// -// Pangle = A*( alpha*PconeA+(1.0-alpha)*PconeB ); -// } -// else{ -// -// // pdfs outside cone -// Pmu = 1.4944765*sin( (Angle0+deltaAngle)*(TMath::Pi()/180.0) ) -// *( 1.0/(1.0+alphaMu*(lambdaMuShort/lambdaMuLong)) )*( 1.0/(1.0+(deltaAngle*deltaAngle)/(lambdaMuShort*lambdaMuShort)) -// + alphaMu*(lambdaMuShort/lambdaMuLong)/(1.0+(deltaAngle*deltaAngle)/(lambdaMuLong*lambdaMuLong)) ); -// -// Pel = 1.4944765*sin( (Angle0+deltaAngle)*(TMath::Pi()/180.0) ) -// *( 1.0/(1.0+alphaEl*(lambdaElShort/lambdaElLong)) )*( 1.0/(1.0+(deltaAngle*deltaAngle)/(lambdaElShort*lambdaElShort)) -// + alphaEl*(lambdaElShort/lambdaElLong)/(1.0+(deltaAngle*deltaAngle)/(lambdaElLong*lambdaElLong)) ); -// -// Pangle = A*( beta*Pmu+(1.0-beta)*Pel ); -// } -// -// // overall probability -// // =================== -// P = Pcharge*Pangle; -// -// chi2 += -2.0*log(P); -// ndof += 1.0; -// } -// } -// -// // calculate figure of merit -// // ========================= -// if( ndof>0.0 ){ -// fom = fBaseFOM - 5.0*chi2/ndof; -// angle = beta*43.0 + (1.0-beta)*49.0; -// } -// -// // return figure of merit -// // ====================== -// coneAngle = angle; -// coneFOM = fom; -// -// return; -//} - +#include "FoMCalculator.h" + +//Constructor +FoMCalculator::FoMCalculator() { + fVtxGeo = 0; + fBaseFOM = 100.0; + + // fitting parameters ported from vertexFinder (FIXME) + fTimeFitWeight = 0.50; // nominal time weight + fConeFitWeight = 0.50; // nominal cone weight + + // default Mean time calculator type + fMeanTimeCalculatorType = 0; +} + +//Destructor +FoMCalculator::~FoMCalculator() { + fVtxGeo = 0; +} + + +void FoMCalculator::LoadVertexGeometry(VertexGeometry* vtxgeo) { + this->fVtxGeo = vtxgeo; +} + + +void FoMCalculator::TimePropertiesLnL(double vtxTime, double& vtxFOM) +{ + // internal variables + // ================== + double weight = 0.0; + double delta = 0.0; // time residual of each hit + double sigma = 0.0; // time resolution of each hit + double A = 0.0; // normalisation of first Gaussian + int type; // Digit type (LAPPD or PMT) + + double Preal = 0.0; // probability of real hit + double P = 0.0; // probability of hit + + double chi2 = 0.0; // log-likelihood: chi2 = -2.0*log(L) + double ndof = 0.0; // total number of hits + double fom = -9999.; // figure of merit + + // add noise to model + // ================== + double Pnoise; + + //FIXME: We need an implementation of noise models for PMTs and LAPPDs + //double nFilterDigits = this->fVtxGeo->GetNFilterDigits(); + //double fTimeFitNoiseRate = 0.02; // hits/ns [0.40 for electrons, 0.02 for muons] + //Pnoise = fTimeFitNoiseRate/nFilterDigits; + + // loop over digits + // ================ + for( int idigit=0; idigitfVtxGeo->GetNDigits(); idigit++ ){ + int detType = this->fVtxGeo->GetDigitType(idigit); + delta = this->fVtxGeo->GetDelta(idigit) - vtxTime; + sigma = this->fVtxGeo->GetDeltaSigma(idigit); + type = this->fVtxGeo->GetDigitType(idigit); + if (type == RecoDigit::PMT8inch){ //PMT8Inch + sigma = 1.5*sigma; + Pnoise = 1e-8; //FIXME; Need implementation of noise model + } + if (type == RecoDigit::lappd_v0){ //lappd + sigma = 1.2*sigma; + Pnoise = 1e-8; //FIXME; Need implementation of noise model + } + A = 1.0 / ( 2.0*sigma*sqrt(0.5*TMath::Pi()) ); //normalisation constant + Preal = A*exp(-(delta*delta)/(2.0*sigma*sigma)); + P = (1.0-Pnoise)*Preal + Pnoise; + chi2 += -2.0*log(P); + ndof += 1.0; + } + + // calculate figure of merit + if( ndof>0.0 ){ + fom = fBaseFOM - 5.0*chi2/ndof; + } + + // return figure of merit + // ====================== + vtxFOM = fom; + return; +} + + + + +void FoMCalculator::ConePropertiesFoM(double coneEdge, double& coneFOM) +{ + // calculate figure of merit + // ========================= + double coneEdgeLow = 21.0; // cone edge (low side) + double coneEdgeHigh = 3.0; // cone edge (high side) [muons: 3.0, electrons: 7.0] + double deltaAngle = 0.0; + double digitCharge = 0.0; + double coneCharge = 0.0; + double allCharge = 0.0; + double outerCone = -99.9; + int outhits = 0; + int inhits = 0; + + double fom = -9999.; + + for( int idigit=0; idigitfVtxGeo->GetNDigits(); idigit++ ){ + if( this->fVtxGeo->IsFiltered(idigit) && this->fVtxGeo->GetDigitType(idigit) == RecoDigit::PMT8inch){ + deltaAngle = this->fVtxGeo->GetAngle(idigit) - coneEdge; + digitCharge = this->fVtxGeo->GetDigitQ(idigit); + + if( deltaAngle<=0.0 ){ + coneCharge += digitCharge*( 0.75 + 0.25/( 1.0 + (deltaAngle*deltaAngle)/(coneEdgeLow*coneEdgeLow) ) ); + inhits++; + //if (deltaAngle > outerCone) outerCone = deltaAngle; + } + else{ + coneCharge += digitCharge*( 0.00 + 1.00/( 1.0 + (deltaAngle*deltaAngle)/(coneEdgeHigh*coneEdgeHigh) ) ); + outhits++; + //outerCone = 0; + } + + allCharge += digitCharge; + //outerCone = -outhits/inhits; + } + } + + if( allCharge>0.0 ){ + if( outerCone>-42 ){ + fom = fBaseFOM*coneCharge/allCharge/*exp(outerCone)*/; + }else{ + fom = fBaseFOM*coneCharge/allCharge; + } + } + + // return figure of merit + // ====================== + coneFOM = fom; + return; +} + +void FoMCalculator::ConePropertiesLnL(double vtxX, double vtxY, double vtxZ, double dirX, double dirY, double dirZ, double coneEdge, double& chi2, TH1D angularDist, double& phimax, double& phimin) { + double coneEdgeLow = 21.0; // cone edge (low side) + double coneEdgeHigh = 3.0; // cone edge (high side) [muons: 3.0, electrons: 7.0] + double deltaAngle = 0.0; + double digitCharge = 0.0; + double digitPE = 0.0; + double coneCharge = 0.0; + double allCharge = 0.0; + double outerCone = -99.9; + double coef = angularDist.Integral(); //1000; + chi2 = 0; + cout << "ConePropertiesLnL Position: (" << vtxX << ", " << vtxY << ", " << vtxZ << ")" << endl; + cout << "And Direction: (" << dirX << ", " << dirY << ", " << dirZ << ")" << endl; + + double digitX, digitY, digitZ; + double dx, dy, dz, ds; + double px, py, pz; + double cosphi, phi, phideg; + phimax = 0; + phimin = 10; + double allPE = 0; + int refbin; + double weight; + double P; + + for (int idigit = 0; idigit < this->fVtxGeo->GetNDigits(); idigit++) { + if (this->fVtxGeo->IsFiltered(idigit) && this->fVtxGeo->GetDigitType(idigit) == RecoDigit::PMT8inch) { + digitCharge = this->fVtxGeo->GetDigitQ(idigit); + allCharge += digitCharge; + } + } + + for (int idigit = 0; idigit < this->fVtxGeo->GetNDigits(); idigit++) { + if (this->fVtxGeo->IsFiltered(idigit) && this->fVtxGeo->GetDigitType(idigit) == RecoDigit::PMT8inch) { + deltaAngle = this->fVtxGeo->GetAngle(idigit) - coneEdge; + digitCharge = this->fVtxGeo->GetDigitQ(idigit); + //digitPE = this->fVtxGeo->GetDigitPE(idigit); + digitX = fVtxGeo->GetDigitX(idigit); + digitY = fVtxGeo->GetDigitY(idigit); + digitZ = fVtxGeo->GetDigitZ(idigit); + dx = digitX - vtxX; + dy = digitY - vtxY; + dz = digitZ - vtxZ; + std::cout << "dx, dy, dz: " << dx << ", " << dy << ", " << dz << endl; + ds = pow(dx * dx + dy * dy + dz * dz, 0.5); + std::cout << "ds: " << ds << endl; + px = dx / ds; + py = dy / ds; + pz = dz / ds; + std::cout << "px, py, pz: " << px << ", " << py << ", " << pz << endl; + std::cout << "dirX, dirY, DirZ: " << dirX << ", " << dirY << ", " << dirZ << endl; + + cosphi = 1.0; + phi = 0.0; + //cout << "angle direction: " << dx << " " << dy << " " << dz << " = " << ds << endl; + cosphi = px * dirX + py * dirY + pz * dirZ; + //cout << "cosphi: " << cosphi << endl; + phi = acos(cosphi); + + if (phi > phimax) phimax = phi; + if (phi < phimin) phimin = phi; + + phideg = phi / (TMath::Pi() / 180); + std::cout << "phi, phideg: " << phi << ", " << phideg << endl; + std::cout << "vs. Zenith: " << fVtxGeo->GetZenith(idigit) << endl; + refbin = angularDist.FindBin(phideg); + weight = angularDist.GetBinContent(refbin)/coef; + P = digitCharge / allCharge; + //cout << "conefomlnl P: " << P << ", weight: " << weight << endl; + chi2 += pow(P - weight, 2)/weight; + + //outerCone = -outhits/inhits; + } + } + //chi2 = (100 - chi2) * exp(-pow(pow(0.7330382, 2) - pow(phimax - phimin, 2), 2) / pow(0.7330382, 2)); +} + + + + +//Given the position of the point vertex (x, y, z) and n digits, calculate the mean expected vertex time +double FoMCalculator::FindSimpleTimeProperties(double myConeEdge) { + double meanTime = 0.0; + // weighted average + if(fMeanTimeCalculatorType == 0) { + // calculate mean and rms of hits inside cone + // ========================================== + double Swx = 0.0; + double Sw = 0.0; + + double delta = 0.0; + double sigma = 0.0; + double weight = 0.0; + double deweight = 0.0; + double deltaAngle = 0.0; + + double myConeEdgeSigma = 7.0; // [degrees] + + for( int idigit=0; idigitfVtxGeo->GetNDigits(); idigit++ ){ + int detType = this->fVtxGeo->GetDigitType(idigit); + if( this->fVtxGeo->IsFiltered(idigit) ){ + delta = this->fVtxGeo->GetDelta(idigit); + sigma = this->fVtxGeo->GetDeltaSigma(idigit); + weight = 1.0/(sigma*sigma); + // profile in angle + deltaAngle = this->fVtxGeo->GetAngle(idigit) - myConeEdge; + // deweight hits outside cone + if( deltaAngle<=0.0 ){ + deweight = 1.0; + } + else{ + deweight = 1.0/(1.0+(deltaAngle*deltaAngle)/(myConeEdgeSigma*myConeEdgeSigma)); + } + Swx += deweight*weight*delta; //delta is expected vertex time + Sw += deweight*weight; + } + } + if( Sw>0.0 ){ + meanTime = Swx*1.0/Sw; + } + } + + // most probable time + else if(fMeanTimeCalculatorType == 1) { + double sigma = 0.0; + double deltaAngle = 0.0; + double weight = 0.0; + double deweight = 0.0; + double myConeEdgeSigma = 7.0; // [degrees] + vector deltaTime1; + vector deltaTime2; + vector TimeWeight; + + for( int idigit=0; idigitGetNDigits(); idigit++ ){ + if(fVtxGeo->IsFiltered(idigit)){ + deltaTime1.push_back(fVtxGeo->GetDelta(idigit)); + deltaTime2.push_back(fVtxGeo->GetDelta(idigit)); + sigma = fVtxGeo->GetDeltaSigma(idigit); + weight = 1.0/(sigma*sigma); + deltaAngle = fVtxGeo->GetAngle(idigit) - myConeEdge; + if( deltaAngle<=0.0 ){ + deweight = 1.0; + } + else{ + deweight = 1.0/(1.0+(deltaAngle*deltaAngle)/(myConeEdgeSigma*myConeEdgeSigma)); + } + TimeWeight.push_back(deweight*weight); + } + } + int n = deltaTime1.size(); + std::sort(deltaTime1.begin(),deltaTime1.end()); + double timeMin = deltaTime1.at(int((n-1)*0.05)); // 5% of the total entries + double timeMax = deltaTime1.at(int((n-1)*0.90)); // 90% of the total entries + int nbins = int(n/5); + TH1D *hDeltaTime = new TH1D("hDeltaTime", "hDeltaTime", nbins, timeMin, timeMax); + for(int i=0; iFill(deltaTime2.at(i), TimeWeight.at(i)); + hDeltaTime->Fill(deltaTime2.at(i)); + } + meanTime = hDeltaTime->GetBinCenter(hDeltaTime->GetMaximumBin()); + delete hDeltaTime; hDeltaTime = 0; + } + + else std::cout<<"FoMCalculator Error: Wrong type of Mean time calculator! "<fVtxGeo->CalcPointResiduals(vtxX, vtxY, vtxZ, 0.0, 0.0, 0.0, 0.0); //calculate expected point vertex time for each digit + + // calculate figure of merit + // ========================= + this->TimePropertiesLnL(vtxTime, vtxFOM); + + // calculate overall figure of merit + // ================================= + fom = vtxFOM; + // truncate + if( fom<-9999. ) fom = -9999.; + + return; +} + + + + + +void FoMCalculator::PointDirectionChi2(double vtxX, double vtxY, + double vtxZ, double dirX, double dirY, double dirZ, + double coneAngle, double& fom) +{ + // figure of merit + // =============== + double coneFOM = -9999.; + + // calculate residuals + // =================== + this->fVtxGeo->CalcPointResiduals(vtxX, vtxY, vtxZ, 0.0, + dirX, dirY, dirZ); //load expected vertex time for each digit + + // calculate figure of merit + // ========================= + this->ConePropertiesFoM(coneAngle, coneFOM); + + // calculate overall figure of merit + // ================================= + fom = coneFOM; + + // truncate + if( fom<-9999. ) fom = -9999.; + + return; +} + +void FoMCalculator::PointVertexChi2(double vtxX, double vtxY, double vtxZ, + double dirX, double dirY, double dirZ, + double coneAngle, double vtxTime, double& fom) +{ + // figure of merit + // =============== + double vtxFOM = -9999.; + + // calculate residuals + // =================== + this->fVtxGeo->CalcPointResiduals(vtxX, vtxY, vtxZ, 0.0, + dirX, dirY, dirZ); //calculate expected vertex time for each digit + // calculate figure of merit + // ========================= + double timeFOM = -9999.; + double coneFOM = -9999.; + + this->ConePropertiesFoM(coneAngle,coneFOM); + this->TimePropertiesLnL(vtxTime, timeFOM); + + double fTimeFitWeight = this->fTimeFitWeight; + double fConeFitWeight = this->fConeFitWeight; + vtxFOM = (fTimeFitWeight*timeFOM+fConeFitWeight*coneFOM)/(fTimeFitWeight+fConeFitWeight); + + // calculate overall figure of merit + // ================================= + fom = vtxFOM; + // truncate + if( fom<-9999. ) fom = -9999.; + + return; +} + +void FoMCalculator::ExtendedVertexChi2(double vtxX, double vtxY, double vtxZ, double dirX, double dirY, double dirZ, double coneAngle, double vtxTime, double& fom) +{ + // figure of merit + // =============== + double vtxFOM = -9999.; + double timeFOM = -9999.; + double coneFOM = -9999.; + + // calculate residuals + // =================== + this->fVtxGeo->CalcExtendedResiduals(vtxX,vtxY,vtxZ,0.0,dirX,dirY,dirZ); + + // calculate figure of merit + // ========================= + + this->ConePropertiesFoM(coneAngle,coneFOM); + this->TimePropertiesLnL(vtxTime, timeFOM); + + double fTimeFitWeight = this->fTimeFitWeight; + double fConeFitWeight = this->fConeFitWeight; + vtxFOM = (fTimeFitWeight*timeFOM+fConeFitWeight*coneFOM)/(fTimeFitWeight+fConeFitWeight); + + // calculate overall figure of merit + // ================================= + fom = vtxFOM; + + // truncate + if( fom<-9999. ) fom = -9999.; + + return; +} + +void FoMCalculator::ExtendedVertexChi2(double vtxX, double vtxY, double vtxZ, double dirX, double dirY, double dirZ, double coneAngle, double vtxTime, double& fom, TH1D pdf) +{ + // figure of merit + // =============== + double vtxFOM = -9999.; + double timeFOM = -9999.; + double coneFOM = -9999.; + double phimax, phimin; + + // calculate residuals + // =================== + this->fVtxGeo->CalcExtendedResiduals(vtxX, vtxY, vtxZ, 0.0, dirX, dirY, dirZ); + + // calculate figure of merit + // ========================= + + this->ConePropertiesLnL(vtxX, vtxY, vtxZ, dirX, dirY, dirZ, coneAngle, coneFOM, pdf, phimax, phimin); + this->TimePropertiesLnL(vtxTime, timeFOM); + + double fTimeFitWeight = this->fTimeFitWeight; + double fConeFitWeight = this->fConeFitWeight; + vtxFOM = (fTimeFitWeight*timeFOM + fConeFitWeight * coneFOM) / (fTimeFitWeight + fConeFitWeight); + + // calculate overall figure of merit + // ================================= + fom = vtxFOM; + + // truncate + if (fom < -9999.) fom = -9999.; + + return; +} + + + +//KEPT FOR HISTORY, BUT FITTER IS CURRENTLY NOT WORKING +//void FoMCalculator::CorrectedVertexChi2(double vtxX, double vtxY, double vtxZ, double dirX, double dirY, double dirZ, double& vtxAngle, double& vtxTime, double& fom) +//{ +// // figure of merit +// // =============== +// double vtxFOM = 0.0; +// double timeFOM = 0.0; +// double coneFOM = 0.0; +// double penaltyFOM = 0.0; +// double fixPositionFOM = 0.0; +// double fixDirectionFOM = 0.0; +// +// // calculate residuals +// // =================== +// this->fVtxGeo->CalcExtendedResiduals(vtxX,vtxY,vtxZ,0.0,dirX,dirY,dirZ); +// +// // calculate figure of merit +// // ========================= +// +// this->ConePropertiesFoM(vtxAngle,coneFOM); +// fom = coneFOM; +// +// // truncate +// if( fom<-999.999*fBaseFOM ) fom = -999.999*fBaseFOM; +// +// return; +//} +// +//void FoMCalculator::ConePropertiesLnL(double coneParam0, double coneParam1, double coneParam2, double& coneAngle, double& coneFOM) +//{ +// +// // nuisance parameters +// // =================== +// double alpha = coneParam0; // track length parameter = 0.25 +// double alpha0 = coneParam1; // track length parameter = 0.5 +// double beta = coneParam2; // particle ID: 0.0[electron]->1.0[muon] = 0.75 +// +// // internal variables +// // ================== +// double deltaAngle = 0.0; // +// double sigmaAngle = 7.0; //Cherenkov Angle resolution +// double Angle0 = Parameters::CherenkovAngle(); //Cherenkov Angle: 42 degree +// double deltaAngle0 = Angle0*alpha0; //? +// +// double digitQ = 0.0; +// double sigmaQmin = 1.0; +// double sigmaQmax = 10.0; +// double sigmaQ = 0.0; +// +// double A = 0.0; +// +// double PconeA = 0.0; +// double PconeB = 0.0; +// double Pmu = 0.0; +// double Pel = 0.0; +// +// double Pcharge = 0.0; +// double Pangle = 0.0; +// double P = 0.0; +// +// double chi2 = 0.0; +// double ndof = 0.0; +// +// double angle = 46.0; +// double fom = 0.0; +// +// // hard-coded parameters: 200 kton (100 kton) +// // ========================================== +// double lambdaMuShort = 0.5; // 0.5; +// double lambdaMuLong = 5.0; // 15.0; +// double alphaMu = 1.0; // 4.5; +// +// double lambdaElShort = 1.0; // 2.5; +// double lambdaElLong = 7.5; // 15.0; +// double alphaEl = 6.0; // 3.5; +// +// // numerical integrals +// // =================== +// fSconeA = 21.9938; +// fSconeB = 0.0000; +// +// // Number of P.E. angular distribution +// // inside cone +// int nbinsInside = 420; //divide the angle range by 420 bins +// for( int n=0; nfVtxGeo->GetNDigits(); idigit++ ){ +// +// if( this->fVtxGeo->IsFiltered(idigit) ){ +// digitQ = this->fVtxGeo->GetDigitQ(idigit); +// deltaAngle = this->fVtxGeo->GetAngle(idigit) - Angle0; +// +// // pulse height distribution +// // ========================= +// if( deltaAngle<=0 ){ //inside Cone +// sigmaQ = sigmaQmax; +// } +// else{ //outside Cone +// sigmaQ = sigmaQmin + (sigmaQmax-sigmaQmin)/(1.0+(deltaAngle*deltaAngle)/(sigmaAngle*sigmaAngle)); +// } +// +// A = 1.0/(log(2.0)+0.5*TMath::Pi()*sigmaQ); +// +// if( digitQ<1.0 ){ +// Pcharge = 2.0*A*digitQ/(1.0+digitQ*digitQ); +// } +// else{ +// Pcharge = A/(1.0+(digitQ-1.0)*(digitQ-1.0)/(sigmaQ*sigmaQ)); +// } +// +// // angular distribution +// // ==================== +// A = 1.0/( alpha*fSconeA + (1.0-alpha)*fSconeB +// + beta*fSmu + (1.0-beta)*fSel ); // numerical integrals +// +// if( deltaAngle<=0 ){ +// +// // pdfs inside cone: +// PconeA = 1.4944765*sin( (Angle0+deltaAngle)*(TMath::Pi()/180.0) ); +// PconeB = 1.4944765*sin( (Angle0+deltaAngle)*(TMath::Pi()/180.0) ) +// *( 1.0/(1.0+(deltaAngle*deltaAngle)/(deltaAngle0*deltaAngle0)) ); +// +// Pangle = A*( alpha*PconeA+(1.0-alpha)*PconeB ); +// } +// else{ +// +// // pdfs outside cone +// Pmu = 1.4944765*sin( (Angle0+deltaAngle)*(TMath::Pi()/180.0) ) +// *( 1.0/(1.0+alphaMu*(lambdaMuShort/lambdaMuLong)) )*( 1.0/(1.0+(deltaAngle*deltaAngle)/(lambdaMuShort*lambdaMuShort)) +// + alphaMu*(lambdaMuShort/lambdaMuLong)/(1.0+(deltaAngle*deltaAngle)/(lambdaMuLong*lambdaMuLong)) ); +// +// Pel = 1.4944765*sin( (Angle0+deltaAngle)*(TMath::Pi()/180.0) ) +// *( 1.0/(1.0+alphaEl*(lambdaElShort/lambdaElLong)) )*( 1.0/(1.0+(deltaAngle*deltaAngle)/(lambdaElShort*lambdaElShort)) +// + alphaEl*(lambdaElShort/lambdaElLong)/(1.0+(deltaAngle*deltaAngle)/(lambdaElLong*lambdaElLong)) ); +// +// Pangle = A*( beta*Pmu+(1.0-beta)*Pel ); +// } +// +// // overall probability +// // =================== +// P = Pcharge*Pangle; +// +// chi2 += -2.0*log(P); +// ndof += 1.0; +// } +// } +// +// // calculate figure of merit +// // ========================= +// if( ndof>0.0 ){ +// fom = fBaseFOM - 5.0*chi2/ndof; +// angle = beta*43.0 + (1.0-beta)*49.0; +// } +// +// // return figure of merit +// // ====================== +// coneAngle = angle; +// coneFOM = fom; +// +// return; +//} diff --git a/DataModel/FoMCalculator.h b/DataModel/FoMCalculator.h old mode 100644 new mode 100755 index a9172ded3..fdd51bdc2 --- a/DataModel/FoMCalculator.h +++ b/DataModel/FoMCalculator.h @@ -1,5 +1,5 @@ /************************************************************************* - > File Name: MinuitOptimizer.hh + > File Name: FoMCalculator.hh > Author: Jingbo Wang, Teal Pershing > mail: jiowang@ucdavis.edu, tjpershing@ucdavis.edu > Created Time: MAY 07, 2019 @@ -50,6 +50,7 @@ class FoMCalculator { double FindSimpleTimeProperties(double myConeEdge); void TimePropertiesLnL(double vtxTime, double& vtxFom); void ConePropertiesFoM(double coneEdge, double& chi2); + void ConePropertiesLnL(double vtxX, double vtxY, double VtxZ, double dirX, double dirY, double dirZ, double coneEdge, double& chi2, TH1D angularDist, double& phimax, double& phimin); void PointPositionChi2(double vtxX, double vtxY, double vtxZ, double vtxTime, double& fom); void PointDirectionChi2(double vtxX, double vtxY, double vtxZ, double dirX, double dirY, double dirZ, double coneAngle, double& fom); void PointVertexChi2(double vtxX, double vtxY, double vtxZ, @@ -58,6 +59,9 @@ class FoMCalculator { void ExtendedVertexChi2(double vtxX, double vtxY, double vtxZ, double dirX, double dirY, double dirZ, double coneAngle, double vtxTime, double& fom); + void ExtendedVertexChi2(double vtxX, double vtxY, double vtxZ, + double dirX, double dirY, double dirZ, + double coneAngle, double vtxTime, double& fom, TH1D pdf); // void ConePropertiesLnL(double coneParam0, double coneParam1, double coneParam2, double& coneAngle, double& coneFOM); // void CorrectedVertexChi2(double vtxX, double vtxY, double vtxZ, // double dirX, double dirY, double dirZ, diff --git a/DataModel/MinuitOptimizer.cpp b/DataModel/MinuitOptimizer.cpp old mode 100644 new mode 100755 index bd9fd7097..34a9911dc --- a/DataModel/MinuitOptimizer.cpp +++ b/DataModel/MinuitOptimizer.cpp @@ -904,6 +904,155 @@ void MinuitOptimizer::FitExtendedVertexWithMinuit() { return; } +void MinuitOptimizer::FitExtendedVertexWithMinuit(TH1D pdf) { + // seed vertex + // =========== + bool foundSeed = (fSeedVtx->FoundVertex() + && fSeedVtx->FoundDirection()); + + + double seedTime = fSeedVtx->GetTime(); + double seedX = fSeedVtx->GetPosition().X(); + double seedY = fSeedVtx->GetPosition().Y(); + double seedZ = fSeedVtx->GetPosition().Z(); + + double seedDirX = fSeedVtx->GetDirection().X(); + double seedDirY = fSeedVtx->GetDirection().Y(); + double seedDirZ = fSeedVtx->GetDirection().Z(); + + double seedTheta = acos(seedDirZ); + double seedPhi = 0.0; + + //modified by JW + if (seedDirX > 0.0) { + seedPhi = atan(seedDirY / seedDirX); + } + if (seedDirX < 0.0) { + seedPhi = atan(seedDirY / seedDirX); + if (seedDirY > 0.0) seedPhi += TMath::Pi(); + if (seedDirY <= 0.0) seedPhi -= TMath::Pi(); + } + if (seedDirX == 0.0) { + if (seedDirY > 0.0) seedPhi = 0.5 * TMath::Pi(); + else if (seedDirY < 0.0) seedPhi = -0.5 * TMath::Pi(); + else seedPhi = 0.0; + } + // current status + // ============== + int status = fSeedVtx->GetStatus(); + + // reset counter + // ============= + extended_vertex_reset_itr(); + + // abort if necessary + // ================== + if (foundSeed == 0) { + if (fPrintLevel >= 0) { + std::cout << " extended vertex fit failed to find input vertex " << std::endl; + } + status |= RecoVertex::kFailExtendedVertex; + fFittedVtx->SetStatus(status); + return; + } + + // run Minuit + // ========== + // six-parameter fit to vertex position, time and direction + + int err = 0; + int flag = 0; + + double fitXpos = 0.0; + double fitYpos = 0.0; + double fitZpos = 0.0; + double fitTheta = 0.0; + double fitPhi = 0.0; + double fitTime = 0.0; + + double fitTimeErr = 0.0; + double fitXposErr = 0.0; + double fitYposErr = 0.0; + double fitZposErr = 0.0; + double fitThetaErr = 0.0; + double fitPhiErr = 0.0; + + double* arglist = new double[10]; + arglist[0] = 2; // 1: standard minimization + // 2: try to improve minimum + +// re-initialize everything... + fMinuitExtendedVertex->mncler(); + fMinuitExtendedVertex->SetFCN(extended_vertex_chi2); + fMinuitExtendedVertex->mnset(); + fMinuitExtendedVertex->mnexcm("SET STR", arglist, 1, err); + fMinuitExtendedVertex->mnparm(0, "x", seedX, 1.0, fXmin, fXmax, err); + fMinuitExtendedVertex->mnparm(1, "y", seedY, 1.0, fYmin, fYmax, err); + fMinuitExtendedVertex->mnparm(2, "z", seedZ, 5.0, fZmin, fZmax, err); + fMinuitExtendedVertex->mnparm(3, "theta", seedTheta, 0.125 * TMath::Pi(), -1.0 * TMath::Pi(), 2.0 * TMath::Pi(), err); + fMinuitExtendedVertex->mnparm(4, "phi", seedPhi, 0.125 * TMath::Pi(), -2.0 * TMath::Pi(), 2.0 * TMath::Pi(), err); + fMinuitExtendedVertex->mnparm(5, "vtxTime", seedTime, 1.0, fTmin, fTmax, err); //....TX + + flag = fMinuitExtendedVertex->Migrad(); + + fMinuitExtendedVertex->GetParameter(0, fitXpos, fitXposErr); + fMinuitExtendedVertex->GetParameter(1, fitYpos, fitYposErr); + fMinuitExtendedVertex->GetParameter(2, fitZpos, fitZposErr); + fMinuitExtendedVertex->GetParameter(3, fitTheta, fitThetaErr); + fMinuitExtendedVertex->GetParameter(4, fitPhi, fitPhiErr); + fMinuitExtendedVertex->GetParameter(5, fitTime, fitTimeErr); + + delete[] arglist; + + //correct angles, JW + if (fitTheta < 0.0) fitTheta = -1.0 * fitTheta; + if (fitTheta > TMath::Pi()) fitTheta = 2.0 * TMath::Pi() - fitTheta; + if (fitPhi < -1.0 * TMath::Pi()) fitPhi += 2.0 * TMath::Pi(); + if (fitPhi > TMath::Pi()) fitPhi -= 2.0 * TMath::Pi(); + + // sort results + // ============ + fVtxX = fitXpos; + fVtxY = fitYpos; + fVtxZ = fitZpos; + fVtxTime = fitTime; + fDirX = sin(fitTheta) * cos(fitPhi); + fDirY = sin(fitTheta) * sin(fitPhi); + fDirZ = cos(fitTheta); + + fVtxFOM = -9999.0; + + fPass = 0; // flag = 0: normal termination + if (flag == 0) fPass = 1; // anything else: abnormal termination + + fItr = extended_vertex_iterations(); + + // fit complete; calculate fit results + // ================ + fgFoMCalculator->ExtendedVertexChi2(fVtxX, fVtxY, fVtxZ, + fDirX, fDirY, fDirZ, + fConeAngle, fVtxTime, fVtxFOM, pdf); + + // set vertex and direction + // ======================== + fFittedVtx->SetVertex(fVtxX, fVtxY, fVtxZ, fVtxTime); + fFittedVtx->SetDirection(fDirX, fDirY, fDirZ); + fFittedVtx->SetConeAngle(fConeAngle); + fFittedVtx->SetFOM(fVtxFOM, fItr, fPass); + + // set status + // ========== + bool inside_det = ANNIEGeometry::Instance()->InsideDetector(fVtxX, fVtxY, fVtxZ); + if (!fPass || !inside_det) status |= RecoVertex::kFailExtendedVertex; + fFittedVtx->SetStatus(status); + if (fPrintLevel >= 0) { + if (!fPass) std::cout << " extended vertex fit failed to converge! Error code: " << flag << std::endl; + } + // return vertex + // ============= + return; +} + //KEPT FOR HISTORY, BUT FITTER IS CURRENTLY NOT WORKING //THESE SHOULD BE MOVED TO BEFORE THE CONSTRUCTOR diff --git a/DataModel/MinuitOptimizer.h b/DataModel/MinuitOptimizer.h old mode 100644 new mode 100755 index 938a48412..fed77588a --- a/DataModel/MinuitOptimizer.h +++ b/DataModel/MinuitOptimizer.h @@ -87,6 +87,7 @@ class MinuitOptimizer { void FitPointDirectionWithMinuit(); void FitPointVertexWithMinuit(); void FitExtendedVertexWithMinuit(); + void FitExtendedVertexWithMinuit(TH1D pdf); double GetTime() {return fVtxTime;} double GetFOM() {return fVtxFOM;} diff --git a/UserTools/Factory/Factory.cpp b/UserTools/Factory/Factory.cpp index cd6659bd5..a5b96e603 100644 --- a/UserTools/Factory/Factory.cpp +++ b/UserTools/Factory/Factory.cpp @@ -154,8 +154,10 @@ if (tool=="checkLAPPDStatus") ret=new checkLAPPDStatus; if (tool=="GetLAPPDEvents") ret=new GetLAPPDEvents; if (tool=="LAPPDDataDecoder") ret=new LAPPDDataDecoder; if (tool=="PythonScript") ret=new PythonScript; +if (tool=="MeanTimeCheck") ret=new MeanTimeCheck; if (tool=="ReweightEventsGenie") ret=new ReweightEventsGenie; if (tool=="FilterLAPPDEvents") ret=new FilterLAPPDEvents; +if (tool=="VtxSeedFineGrid") ret=new VtxSeedFineGrid; if (tool=="FilterEvents") ret=new FilterEvents; if (tool=="Stage1DataBuilder") ret=new Stage1DataBuilder; if (tool=="BeamFetcherV2") ret=new BeamFetcherV2; diff --git a/UserTools/LikelihoodFitterCheck/LikelihoodFitterCheck.cpp b/UserTools/LikelihoodFitterCheck/LikelihoodFitterCheck.cpp old mode 100644 new mode 100755 index 91cda74ef..12c518576 --- a/UserTools/LikelihoodFitterCheck/LikelihoodFitterCheck.cpp +++ b/UserTools/LikelihoodFitterCheck/LikelihoodFitterCheck.cpp @@ -14,14 +14,21 @@ bool LikelihoodFitterCheck::Initialise(std::string configfile, DataModel &data){ m_variables.Get("OutputFile", output_filename); m_variables.Get("ifPlot2DFOM", ifPlot2DFOM); m_variables.Get("ShowEvent", fShowEvent); + m_variables.Get("UsePDFFile", fUsePDFFile); + m_variables.Get("PDFFile", pdffile); fOutput_tfile = new TFile(output_filename.c_str(), "recreate"); // Histograms Likelihood2D = new TH2D("Likelihood2D","Figure of merit 2D", 200, -50, 150, 100, -50, 50); + Likelihood2D_pdf = new TH2D("Likelihood2D_pdf", "pdf-based figure of merit 2D", 200, 0, 200, 100, 0, 100); gr_parallel = new TGraph(); gr_parallel->SetTitle("Figure of merit parallel to the track direction"); gr_transverse = new TGraph(); gr_transverse->SetTitle("Figure of merit transverse to the track direction"); + pdf_parallel = new TGraph(); + pdf_parallel->SetTitle("PDF-based Figure of merit parallel to track direction"); + pdf_transverse = new TGraph(); + pdf_transverse->SetTitle("PDF-based figure of merit transverse to track direction"); m_data= &data; //assigning transient data pointer ///////////////////////////////////////////////////////////////// @@ -30,6 +37,7 @@ bool LikelihoodFitterCheck::Initialise(std::string configfile, DataModel &data){ bool LikelihoodFitterCheck::Execute(){ + Log("===========================================================================================",v_debug,verbosity); Log("LikelihoodFitterCheck Tool: Executing",v_debug,verbosity); @@ -73,6 +81,15 @@ bool LikelihoodFitterCheck::Execute(){ Log("LikelihoodFitterCheck Tool: Error retrieving RecoDigits,no digit from the RecoEvent!",v_error,verbosity); return false; } + + if (fUsePDFFile) { + bool pdftest = this->GetPDF(pdf); + if (!pdftest) { + Log("LikelihoodFitterCheck Tool: Error retrieving pdffile; running without!", v_error, verbosity); + fUsePDFFile = 0; + return false; + } + } double recoVtxX, recoVtxY, recoVtxZ, recoVtxT, recoDirX, recoDirY, recoDirZ; double trueVtxX, trueVtxY, trueVtxZ, trueVtxT, trueDirX, trueDirY, trueDirZ; @@ -101,12 +118,13 @@ bool LikelihoodFitterCheck::Execute(){ double dx = dl * trueDirX; double dy = dl * trueDirY; double dz = dl * trueDirZ; - int nbins = 200; - double dlpara[200], dlfom[200]; - for(int j=0;j<200;j++) { - seedX = trueVtxX - 50*dx + j*dx; - seedY = trueVtxY - 50*dy + j*dy; - seedZ = trueVtxZ - 50*dz + j*dz; + int nbins = 400; + double dlpara[400], dlfom[400]; + double minphi, maxphi; + for(int j=0;j<400;j++) { + seedX = trueVtxX - 350*dx + j*dx; + seedY = trueVtxY - 350*dy + j*dy; + seedZ = trueVtxZ - 350*dz + j*dz; seedT = trueVtxT; seedDirX = trueDirX; seedDirY = trueDirY; @@ -117,14 +135,23 @@ bool LikelihoodFitterCheck::Execute(){ Double_t fom = -999.999*100; double timefom = -999.999*100; double conefom = -999.999*100; + double conefomlnl = -999.999 * 100; + Double_t fompdf = -999.999 * 100; myFoMCalculator->TimePropertiesLnL(meantime,timefom); myFoMCalculator->ConePropertiesFoM(ConeAngle,conefom); fom = timefom*0.5+conefom*0.5; cout<<"timeFOM, coneFOM, fom = "<SetPoint(j, dlpara[j], dlfom[j]); + + if (fUsePDFFile) { + myFoMCalculator->ConePropertiesLnL(seedX, seedY, seedZ, seedDirX, seedDirY, seedDirZ, ConeAngle, conefomlnl, pdf, maxphi, minphi); + cout << "conefomlnl: " << conefomlnl << endl; + fompdf = 0.5 * timefom + 0.5 * conefomlnl; + pdf_parallel->SetPoint(j, dlpara[j], fompdf); + } } //transverse direction @@ -148,8 +175,10 @@ bool LikelihoodFitterCheck::Execute(){ int nhits = myvtxgeo->GetNDigits(); double meantime = myFoMCalculator->FindSimpleTimeProperties(ConeAngle); Double_t fom = -999.999*100; + Double_t fompdf = -999.999 * 100; double timefom = -999.999*100; double conefom = -999.999*100; + double conefomlnl = -999.999 * 100; double coneAngle = 42.0; myFoMCalculator->TimePropertiesLnL(meantime,timefom); myFoMCalculator->ConePropertiesFoM(ConeAngle,conefom); @@ -159,10 +188,17 @@ bool LikelihoodFitterCheck::Execute(){ dltrans[j] = - 50*dl + j*dl; dlfom[j] = fom; gr_transverse->SetPoint(j, dlpara[j], dlfom[j]); + if (fUsePDFFile) { + cout << "pdf fom coming\n"; + // myFoMCalculator->ConePropertiesLnL(seedX, seedY, seedZ, trueDirX, trueDirY, trueDirZ, coneAngle, conefomlnl, pdf); + fompdf = timefom * conefomlnl; + pdf_transverse->SetPoint(j, dltrans[j], conefomlnl); + } } if(ifPlot2DFOM) { //2D scan around the true vertex position + cout << "2DPlot starting now" << endl; double dl_para = 1.0, dl_trans = 1.0; double dx_para = dl_para * trueDirX; double dy_para = dl_para * trueDirY; @@ -170,42 +206,76 @@ bool LikelihoodFitterCheck::Execute(){ double dx_trans = dl_trans * v.X(); double dy_trans = dl_trans * v.Y(); double dz_trans = dl_trans * v.Z(); + double phimax, phimin; for(int k=0; k<100; k++) { for(int m=0; m<200; m++) { seedX = trueVtxX - 50*dx_trans + k*dx_trans - 50*dx_para + m*dx_para; seedY = trueVtxY - 50*dy_trans + k*dy_trans - 50*dy_para + m*dy_para; seedZ = trueVtxZ - 50*dz_trans + k*dz_trans - 50*dz_para + m*dz_para; seedT = trueVtxT; - seedDirX = trueDirX; - seedDirY = trueDirY; - seedDirZ = trueDirZ; - myvtxgeo->CalcExtendedResiduals(seedX, seedY, seedZ, seedT, seedDirX, seedDirY, seedDirZ); + seedDirX = cos(m * TMath::Pi() / 100) * sin(k * TMath::Pi() / 100); + seedDirY = sin(m * TMath::Pi() / 100) * sin(k * TMath::Pi() / 100); + seedDirZ = cos(k * TMath::Pi() / 100); + myvtxgeo->CalcExtendedResiduals(trueVtxX, trueVtxY, trueVtxZ, seedT, seedDirX, seedDirY, seedDirZ); int nhits = myvtxgeo->GetNDigits(); double meantime = myFoMCalculator->FindSimpleTimeProperties(ConeAngle); Double_t fom = -999.999*100; + Double_t fompdf = -999.999 * 100; double timefom = -999.999*100; double conefom = -999.999*100; + double conefomlnl = -999.999 * 100; double coneAngle = 42.0; myFoMCalculator->TimePropertiesLnL(meantime,timefom); myFoMCalculator->ConePropertiesFoM(coneAngle,conefom); - fom = timefom*0.5+conefom*0.5; + fom = timefom * 0.5 + conefom * 0.5; //fom = timefom; cout<<"k,m, timeFOM, coneFOM, fom = "<SetBinContent(m, k, fom); + if (fUsePDFFile) { + myFoMCalculator->ConePropertiesLnL(trueVtxX, trueVtxY, trueVtxZ, seedDirX, seedDirY, seedDirZ, coneAngle, conefomlnl, pdf, phimax, phimin); + fompdf = 0.5 * timefom + 0.5 * (conefomlnl); + cout << "coneFOMlnl: " << conefomlnl << endl; + if (k == 50 && m == 50) { + std::cout << "!!!OUTPUT!!! at true:\n"; + } + std::cout<<"conefomlnl, timefom, fompdf: " << conefomlnl << ", " << timefom << ", " << fompdf << endl; + std::cout << "phimax, phimin: " << phimax << ", " << phimin << endl; + Likelihood2D_pdf->SetBinContent(m, k, fompdf); + } } } } + delete myFoMCalculator; return true; } -bool LikelihoodFitterCheck::Finalise(){ - fOutput_tfile->cd(); - gr_parallel->Write(); - gr_transverse->Write(); - fOutput_tfile->Write(); - fOutput_tfile->Close(); - Log("LikelihoodFitterCheck exitting", v_debug,verbosity); - return true; +bool LikelihoodFitterCheck::Finalise() { + fOutput_tfile->cd(); + gr_parallel->Write(); + gr_transverse->Write(); + + if (fUsePDFFile) { + pdf_parallel->Write(); + pdf_transverse->Write(); + } + + fOutput_tfile->Write(); + fOutput_tfile->Close(); + + Log("LikelihoodFitterCheck exitting", v_debug, verbosity); + return true; } + + + bool LikelihoodFitterCheck::GetPDF(TH1D & pdf) { + TFile f1(pdffile.c_str(), "READ"); + if (!f1.IsOpen()) { + Log("VtxExtendedVertexFinder: pdffile does not exist", v_error, verbosity); + return false; + } + pdf = *(TH1D*)f1.Get("zenith"); + cout << "pdf entries: " << pdf.GetEntries() << endl; + return true; + } diff --git a/UserTools/LikelihoodFitterCheck/LikelihoodFitterCheck.h b/UserTools/LikelihoodFitterCheck/LikelihoodFitterCheck.h old mode 100644 new mode 100755 index 22ed724aa..845629ce7 --- a/UserTools/LikelihoodFitterCheck/LikelihoodFitterCheck.h +++ b/UserTools/LikelihoodFitterCheck/LikelihoodFitterCheck.h @@ -7,6 +7,7 @@ #include "Tool.h" #include "TH1.h" #include "TH2D.h" +#include "TFile.h" #include "FoMCalculator.h" #include "VertexGeometry.h" @@ -22,6 +23,7 @@ class LikelihoodFitterCheck: public Tool { bool Initialise(std::string configfile,DataModel &data); bool Execute(); bool Finalise(); + bool GetPDF(TH1D & pdf); private: @@ -50,6 +52,11 @@ class LikelihoodFitterCheck: public Tool { TH2D* Likelihood2D = 0; TGraph *gr_parallel = 0; TGraph *gr_transverse = 0; + + /// \comparison histograms + TH2D* Likelihood2D_pdf = 0; + TGraph* pdf_parallel = 0; + TGraph* pdf_transverse = 0; /// verbosity levels: if 'verbosity' < this level, the message type will be logged. int verbosity=-1; @@ -60,6 +67,9 @@ class LikelihoodFitterCheck: public Tool { std::string logmessage; int get_ok; bool ifPlot2DFOM = false; + std::string pdffile; + bool fUsePDFFile = 0; + TH1D pdf; diff --git a/UserTools/LikelihoodFitterCheck/README.md b/UserTools/LikelihoodFitterCheck/README.md old mode 100644 new mode 100755 diff --git a/UserTools/MeanTimeCheck/MeanTimeCheck.cpp b/UserTools/MeanTimeCheck/MeanTimeCheck.cpp new file mode 100644 index 000000000..65f6a4ef8 --- /dev/null +++ b/UserTools/MeanTimeCheck/MeanTimeCheck.cpp @@ -0,0 +1,185 @@ +#include "MeanTimeCheck.h" + +MeanTimeCheck::MeanTimeCheck():Tool(){} + + +bool MeanTimeCheck::Initialise(std::string configfile, DataModel &data){ + + /////////////////// Useful header /////////////////////// + if(configfile!="") m_variables.Initialise(configfile); // loading config file + //m_variables.Print(); + + std::string output_filename; + m_variables.Get("verbosity", verbosity); + m_variables.Get("OutputFile", output_filename); + m_variables.Get("ShowEvent", ShowEvent); + Output_tfile = new TFile(output_filename.c_str(), "recreate"); + + delta = new TH1D("delta", "delta", 1000, -10, 30); + peakTime = new TH1D("Peak Time", "time", 1000, -10, 30); + basicAverage = new TH1D("basic average", "time", 1000, -10, 30); + weightedPeak = new TH1D("Weighted Peak", "time", 1000, -10, 30); + weightedAverage = new TH1D("Weighted Average", "time", 1000, -10, 30); + + m_data= &data; //assigning transient data pointer + ///////////////////////////////////////////////////////////////// + + return true; +} + + +bool MeanTimeCheck::Execute(){ + +Log("===========================================================================================", v_debug, verbosity); + + Log("MeanTimeCheck Tool: Executing", v_debug, verbosity); + auto* reco_event = m_data->Stores["RecoEvent"]; + if (!reco_event) { + Log("Error: The MeanTimeCheck tool could not find the ANNIEEvent Store", 0, verbosity); + return false; + } + + m_data->Stores.at("ANNIEEvent")->Get("MCEventNum", MCEventNum); + m_data->Stores.at("ANNIEEvent")->Get("MCTriggernum", MCTriggerNum); + m_data->Stores.at("ANNIEEvent")->Get("EventNumber", EventNumber); + + //std::cout< 0 && EventNumber != ShowEvent) return true; + bool EventCutstatus = false; + auto get_evtstatus = m_data->Stores.at("RecoEvent")->Get("EventCutStatus", EventCutstatus); + if (!get_evtstatus) { + Log("Error: The MeanTimeCheck tool could not find the Event selection status", v_error, verbosity); + return false; + } + if (!EventCutstatus) { + Log("Message: This event doesn't pass the event selection. ", v_message, verbosity); + return true; + } + + auto get_vtx = m_data->Stores.at("RecoEvent")->Get("TrueVertex", TrueVtx); + auto get_digit = m_data->Stores.at("RecoEvent")->Get("RecoDigit", DigitList); + Position vtxPos = TrueVtx->GetPosition(); + Direction vtxDir = TrueVtx->GetDirection(); + double vtxT = TrueVtx->GetTime(); + VtxGeo = VertexGeometry::Instance(); + VtxGeo->LoadDigits(DigitList); + VtxGeo->CalcExtendedResiduals(vtxPos.X(), vtxPos.Y(), vtxPos.Z(), vtxT, vtxDir.X(), vtxDir.Y(), vtxDir.Z()); + int nhits = VtxGeo->GetNDigits(); + //std::cout<<"nhits: "<Fill(VtxGeo->GetDelta(n)); + //std::cout<<"Delta: "<GetDelta(n)<GetPeakTime(); + double WeightedAverage = this->GetWeightedAverage(); + double WeightedPeak = this-> GetWeightedPeak(); +std::cout<<"Check 2"<Fill(PeakTime); + weightedAverage->Fill(WeightedAverage); + weightedPeak->Fill(WeightedPeak); +std::cout<<"Check 3"<Stores.at("RecoEvent")->Set("meanTime", PeakTime); + + return true; + +} + + +bool MeanTimeCheck::Finalise(){ + + Output_tfile->cd(); + Output_tfile->Write(); + Output_tfile->Close(); + Log("MeanTimeCheck exitting", v_debug, verbosity); + + return true; +} + +double MeanTimeCheck::GetPeakTime() { + return delta->GetBinCenter(delta->GetMaximumBin()); +} + +double MeanTimeCheck::GetWeightedAverage() { + double Swx = 0.0; + double Sw = 0.0; + + double delta = 0.0; + double sigma = 0.0; + double weight = 0.0; + double deweight = 0.0; + double deltaAngle = 0.0; + double meanTime = 0.0; + + double myConeEdgeSigma = 7.0; // [degrees] + + for (int idigit = 0; idigit < this->VtxGeo->GetNDigits(); idigit++) { + int detType = this->VtxGeo->GetDigitType(idigit); + if (this->VtxGeo->IsFiltered(idigit)) { + delta = this->VtxGeo->GetDelta(idigit); + sigma = this->VtxGeo->GetDeltaSigma(idigit); + weight = 1.0 / (sigma*sigma); + // profile in angle + deltaAngle = this->VtxGeo->GetAngle(idigit) - Parameters::CherenkovAngle(); + // deweight hits outside cone + if (deltaAngle <= 0.0) { + deweight = 1.0; + } + else { + deweight = 1.0 / (1.0 + (deltaAngle*deltaAngle) / (myConeEdgeSigma*myConeEdgeSigma)); + } + Swx += deweight * weight*delta; //delta is expected vertex time + Sw += deweight * weight; + } + } + if (Sw > 0.0) { + meanTime = Swx * 1.0 / Sw; + } + +return meanTime; +} + +double MeanTimeCheck::GetWeightedPeak() { + double sigma = 0.0; + double deltaAngle = 0.0; + double weight = 0.0; + double deweight = 0.0; + double myConeEdgeSigma = 7.0; // [degrees] + double meanTime = 0.0; + vector deltaTime1; + vector deltaTime2; + vector TimeWeight; + + for (int idigit = 0; idigit < this->VtxGeo->GetNDigits(); idigit++) { + if (this->VtxGeo->IsFiltered(idigit)) { + deltaTime1.push_back(this->VtxGeo->GetDelta(idigit)); + deltaTime2.push_back(this->VtxGeo->GetDelta(idigit)); + sigma = this->VtxGeo->GetDeltaSigma(idigit); + weight = 1.0 / (sigma*sigma); + deltaAngle = this->VtxGeo->GetAngle(idigit) - Parameters::CherenkovAngle(); + if (deltaAngle <= 0.0) { + deweight = 1.0; + } + else { + deweight = 1.0 / (1.0 + (deltaAngle*deltaAngle) / (myConeEdgeSigma*myConeEdgeSigma)); + } + TimeWeight.push_back(deweight*weight); + } + } + int n = deltaTime1.size(); + std::sort(deltaTime1.begin(), deltaTime1.end()); + double timeMin = deltaTime1.at(int((n - 1)*0.05)); // 5% of the total entries + double timeMax = deltaTime1.at(int((n - 1)*0.90)); // 90% of the total entries + int nbins = int(n / 5); + TH1D *hDeltaTime = new TH1D("hDeltaTime", "hDeltaTime", nbins, timeMin, timeMax); + for (int i = 0; i < n; i++) { + hDeltaTime->Fill(deltaTime2.at(i), TimeWeight.at(i)); + hDeltaTime->Fill(deltaTime2.at(i)); + } + meanTime = hDeltaTime->GetBinCenter(hDeltaTime->GetMaximumBin()); + delete hDeltaTime; hDeltaTime = 0; + + return meanTime; +} diff --git a/UserTools/MeanTimeCheck/MeanTimeCheck.h b/UserTools/MeanTimeCheck/MeanTimeCheck.h new file mode 100644 index 000000000..3231ad881 --- /dev/null +++ b/UserTools/MeanTimeCheck/MeanTimeCheck.h @@ -0,0 +1,70 @@ +#ifndef MeanTimeCheck_H +#define MeanTimeCheck_H + +#include +#include +#include "TROOT.h" +#include "TChain.h" +#include "TFile.h" + +#include "Tool.h" +#include "TTree.h" +#include "TH2D.h" +#include "Parameters.h" +#include "VertexGeometry.h" +#include "Detector.h" + + +/** + * \class MeanTimeCheck + * + * This is a blank template for a Tool used by the script to generate a new custom tool. Please fill out the description and author information. +* +* $Author: B.Richards $ +* $Date: 2019/05/28 10:44:00 $ +* Contact: b.richards@qmul.ac.uk +*/ +class MeanTimeCheck: public Tool { + + + public: + + MeanTimeCheck(); ///< Simple constructor + bool Initialise(std::string configfile,DataModel &data); ///< Initialise Function for setting up Tool resources. @param configfile The path and name of the dynamic configuration file to read in. @param data A reference to the transient data class used to pass information between Tools. + bool Execute(); ///< Execute function used to perform Tool purpose. + bool Finalise(); ///< Finalise function used to clean up resources. + + + private: + +double BasicAverage(); +std::vector* DigitList = 0; +RecoVertex* TrueVtx = 0; +double GetPeakTime(); +double GetWeightedAverage(); +double GetWeightedPeak(); +TFile* Output_tfile = nullptr; +TTree* fVertexGeometry = nullptr; +uint64_t MCEventNum; +uint16_t MCTriggerNum; +uint32_t EventNumber; +TH1D *delta; +TH1D *peakTime; +TH1D *basicAverage; +TH1D *weightedAverage; +TH1D *weightedPeak; +VertexGeometry* VtxGeo; + +int verbosity = -1; +int ShowEvent = -1; +int v_error = 0; +int v_warning = 1; +int v_message = 2; +int v_debug = 3; +std::string logmessage; +int get_ok; + +}; + + +#endif diff --git a/UserTools/MeanTimeCheck/README.md b/UserTools/MeanTimeCheck/README.md new file mode 100644 index 000000000..81ae8f285 --- /dev/null +++ b/UserTools/MeanTimeCheck/README.md @@ -0,0 +1,20 @@ +# MeanTimeCheck + +MeanTimeCheck + +## Data + +Describe any data formats MeanTimeCheck creates, destroys, changes, or analyzes. E.G. + +**RawLAPPDData** `map>>` +* Takes this data from the `ANNIEEvent` store and finds the number of peaks + + +## Configuration + +Describe any configuration variables for MeanTimeCheck. + +``` +param1 value1 +param2 value2 +``` diff --git a/UserTools/Unity.h b/UserTools/Unity.h index dc6e501b7..6994cfab4 100644 --- a/UserTools/Unity.h +++ b/UserTools/Unity.h @@ -162,8 +162,10 @@ #include "GetLAPPDEvents.h" #include "LAPPDDataDecoder.h" #include "PythonScript.h" +#include "MeanTimeCheck.h" #include "ReweightEventsGenie.h" #include "FilterLAPPDEvents.h" +#include "VtxSeedFineGrid.h" #include "FilterEvents.h" #include "Stage1DataBuilder.h" #include "BeamFetcherV2.h" diff --git a/UserTools/VertexGeometryCheck/README.md b/UserTools/VertexGeometryCheck/README.md old mode 100644 new mode 100755 diff --git a/UserTools/VertexGeometryCheck/VertexGeometryCheck.cpp b/UserTools/VertexGeometryCheck/VertexGeometryCheck.cpp old mode 100644 new mode 100755 index 9cde9d960..3398f0727 --- a/UserTools/VertexGeometryCheck/VertexGeometryCheck.cpp +++ b/UserTools/VertexGeometryCheck/VertexGeometryCheck.cpp @@ -12,6 +12,8 @@ bool VertexGeometryCheck::Initialise(std::string configfile, DataModel &data){ m_variables.Get("verbosity", verbosity); m_variables.Get("OutputFile", output_filename); m_variables.Get("ShowEvent", fShowEvent); + m_variables.Get("Theta", vertheta); + m_variables.Get("Phi", verphi); fOutput_tfile = new TFile(output_filename.c_str(), "recreate"); // Histograms @@ -106,6 +108,12 @@ bool VertexGeometryCheck::Execute(){ trueDirX = vtxDir.X(); trueDirY = vtxDir.Y(); trueDirZ = vtxDir.Z(); + if (vertheta != -999 && verphi != -999) { + Log("overriding direction", v_debug, verbosity); + trueDirX = cos(vertheta) * sin(verphi); + trueDirY = sin(vertheta) * sin(verphi); + trueDirZ = cos(verphi); + } double ConeAngle = Parameters::CherenkovAngle(); @@ -160,7 +168,7 @@ bool VertexGeometryCheck::Execute(){ if(myvtxgeo->GetDigitType(n)==RecoDigit::PMT8inch) fpmtextendedtres->Fill(myvtxgeo->GetExtendedResidual(n)); fltrack->Fill(myvtxgeo->GetDistTrack(n)); //cm flphoton->Fill(myvtxgeo->GetDistPhoton(n)); //cm - fzenith->Fill(myvtxgeo->GetZenith(n)); // + fzenith->Fill(phideg/*myvtxgeo->GetZenith(n)*/); // fazimuth->Fill(myvtxgeo->GetAzimuth(n)); // fconeangle->Fill(myvtxgeo->GetConeAngle(n)); // fdigitcharge->Fill(myvtxgeo->GetDigitQ(n)); diff --git a/UserTools/VertexGeometryCheck/VertexGeometryCheck.h b/UserTools/VertexGeometryCheck/VertexGeometryCheck.h old mode 100644 new mode 100755 index 414f309b0..41038de58 --- a/UserTools/VertexGeometryCheck/VertexGeometryCheck.h +++ b/UserTools/VertexGeometryCheck/VertexGeometryCheck.h @@ -65,6 +65,7 @@ class VertexGeometryCheck: public Tool { TH1D *flappdtimesmear; TH1D *fpmttimesmear; TH2D *fYvsDigitTheta_all; + double vertheta = -999, verphi = -999; /// verbosity levels: if 'verbosity' < this level, the message type will be logged. diff --git a/UserTools/VtxExtendedVertexFinder/README.md b/UserTools/VtxExtendedVertexFinder/README.md old mode 100644 new mode 100755 diff --git a/UserTools/VtxExtendedVertexFinder/VtxExtendedVertexFinder.cpp b/UserTools/VtxExtendedVertexFinder/VtxExtendedVertexFinder.cpp old mode 100644 new mode 100755 index 27a527059..84b36b685 --- a/UserTools/VtxExtendedVertexFinder/VtxExtendedVertexFinder.cpp +++ b/UserTools/VtxExtendedVertexFinder/VtxExtendedVertexFinder.cpp @@ -22,6 +22,8 @@ bool VtxExtendedVertexFinder::Initialise(std::string configfile, DataModel &data m_variables.Get("verbosity", verbosity); m_variables.Get("FitTimeWindowMin", fTmin); m_variables.Get("FitTimeWindowMax", fTmax); + m_variables.Get("UsePDFFile", fUsePDFFile); + m_variables.Get("PDFFile", pdffile); /// Create extended vertex /// Note that the objects created by "new" must be added to the "RecoEvent" store. @@ -141,7 +143,8 @@ bool VtxExtendedVertexFinder::Finalise(){ } RecoVertex* VtxExtendedVertexFinder::FitExtendedVertex(RecoVertex* myVertex) { - //fit with Minuit + + //fit with Minuit MinuitOptimizer* myOptimizer = new MinuitOptimizer(); myOptimizer->SetPrintLevel(-1); myOptimizer->SetMeanTimeCalculatorType(1); //Type 1: most probable time @@ -177,6 +180,14 @@ RecoVertex* VtxExtendedVertexFinder::FitGridSeeds(std::vector* vSeed RecoVertex* fSimpleVertex = new RecoVertex(); RecoVertex* bestGridVertex = new RecoVertex(); // FIXME: pointer must be deleted by the invoker + if (fUsePDFFile) { + bool pdftest = this->GetPDF(pdf); + if (!pdftest) { + Log("pdffile error; continuing with fom reconstruction", v_error, verbosity); + fUsePDFFile = 0; + } + } + for( unsigned int n=0; n* vSeed fSeedPos = &(vSeedVtxList->at(n)); fSimpleVertex= this->FindSimpleDirection(fSeedPos); myOptimizer->LoadVertex(fSimpleVertex); //Load vertex seed - myOptimizer->FitExtendedVertexWithMinuit(); //scan the point position in 4D space + if (!fUsePDFFile) { + myOptimizer->FitExtendedVertexWithMinuit(); //scan the point position in 4D space + } + else { + myOptimizer->FitExtendedVertexWithMinuit(pdf); + } + vtxFOM = myOptimizer->GetFittedVertex()->GetFOM(); vtxRecoStatus = myOptimizer->GetFittedVertex()->GetStatus(); @@ -203,6 +220,7 @@ RecoVertex* VtxExtendedVertexFinder::FitGridSeeds(std::vector* vSeed std::cout << "best fit reco status: " << bestGridVertex->GetStatus() << std::endl; std::cout << "BestVertex info: " << bestGridVertex->Print() << std::endl; } + return bestGridVertex; } @@ -298,6 +316,16 @@ void VtxExtendedVertexFinder::PushExtendedVertex(RecoVertex* vtx, bool savetodis m_data->Stores.at("RecoEvent")->Set("ExtendedVertex", fExtendedVertex, savetodisk); } +bool VtxExtendedVertexFinder::GetPDF(TH1D& pdf) { + TFile f1(pdffile.c_str(), "READ"); + if (!f1.IsOpen()) { + Log("VtxExtendedVertexFinder: pdffile does not exist", v_error, verbosity); + return false; + } + pdf = *(TH1D*)f1.Get("zenith"); + return true; +} + void VtxExtendedVertexFinder::Reset() { fExtendedVertex->Reset(); } diff --git a/UserTools/VtxExtendedVertexFinder/VtxExtendedVertexFinder.h b/UserTools/VtxExtendedVertexFinder/VtxExtendedVertexFinder.h old mode 100644 new mode 100755 index 92ff7048f..840fd2c58 --- a/UserTools/VtxExtendedVertexFinder/VtxExtendedVertexFinder.h +++ b/UserTools/VtxExtendedVertexFinder/VtxExtendedVertexFinder.h @@ -36,6 +36,12 @@ class VtxExtendedVertexFinder: public Tool { // \brief maximum of fit time window double fTmax; + // \use external file for PDF? If 0, will use equation fit + bool fUsePDFFile = 0; + + // \file containing histogram of PDF of charge-angle distribution + std::string pdffile; + /// \brief RecoVertex* FitExtendedVertex(RecoVertex* myvertex); @@ -44,6 +50,7 @@ class VtxExtendedVertexFinder: public Tool { /// \brief Find a simple direction using weighted sum of digit charges RecoVertex* FindSimpleDirection(RecoVertex* myvertex); + bool GetPDF(TH1D &pdf); /// \brief Reset everything void Reset(); @@ -70,7 +77,8 @@ class VtxExtendedVertexFinder: public Tool { int v_message=2; int v_debug=3; std::string logmessage; - int get_ok; + int get_ok; + TH1D pdf; diff --git a/UserTools/VtxSeedFineGrid/README.md b/UserTools/VtxSeedFineGrid/README.md new file mode 100755 index 000000000..fc5f3d780 --- /dev/null +++ b/UserTools/VtxSeedFineGrid/README.md @@ -0,0 +1,20 @@ +# VtxSeedFineGrid + +VtxSeedFineGrid + +## Data + +Describe any data formats VtxSeedFineGrid creates, destroys, changes, or analyzes. E.G. + +**RawLAPPDData** `map>>` +* Takes this data from the `ANNIEEvent` store and finds the number of peaks + + +## Configuration + +Describe any configuration variables for VtxSeedFineGrid. + +``` +param1 value1 +param2 value2 +``` diff --git a/UserTools/VtxSeedFineGrid/VtxSeedFineGrid.cpp b/UserTools/VtxSeedFineGrid/VtxSeedFineGrid.cpp new file mode 100755 index 000000000..3b69ef17d --- /dev/null +++ b/UserTools/VtxSeedFineGrid/VtxSeedFineGrid.cpp @@ -0,0 +1,412 @@ +#include "VtxSeedFineGrid.h" + +VtxSeedFineGrid::VtxSeedFineGrid():Tool(){} + + +bool VtxSeedFineGrid::Initialise(std::string configfile, DataModel &data){ + + /////////////////// Useful header /////////////////////// + if(configfile!="") m_variables.Initialise(configfile); // loading config file + //m_variables.Print(); + m_variables.Get("verbosity", verbosity); + m_variables.Get("useTrueDir", useTrueDir); + m_variables.Get("useSimpleDir", useSimpleDir); + m_variables.Get("useMRDTrack", useMRDTrack); + m_variables.Get("usePastResolution", usePastResolution); + m_variables.Get("useDirectionGrid", useDirectionGrid); + m_variables.Get("multiGrid", multiGrid); + m_variables.Get("InputFile", InputFile); + + + m_data= &data; //assigning transient data pointer + ///////////////////////////////////////////////////////////////// + vSeedVtxList = nullptr; + + return true; +} + + +bool VtxSeedFineGrid::Execute(){ + Log("VtxSeedFineGrid Tool: Executing", v_debug, verbosity); + + auto* annie_event = m_data->Stores["RecoEvent"]; + if (!annie_event) { + Log("Error: VtxSeedFineGrid tool could not find the RecoEvent Store", v_error, verbosity); + return false; + } + + m_data->Stores.at("RecoEvent")->Get("vSeedVtxList", vSeedVtxList); + + auto get_vtx = m_data->Stores.at("RecoEvent")->Get("TrueVertex", fTrueVertex); ///> Get digits from "RecoEvent" +/* if (!get_vtx) { + Log("LikelihoodFitterCheck Tool: Error retrieving TrueVertex! ", v_error, verbosity); + return false; + }*/ + + // Retrive digits from RecoEvent + auto get_digit = m_data->Stores.at("RecoEvent")->Get("RecoDigit", fDigitList); ///> Get digits from "RecoEvent" + if (!get_digit) { + Log("LikelihoodFitterCheck Tool: Error retrieving RecoDigits,no digit from the RecoEvent!", v_error, verbosity); + return false; + } + auto get_flagsapp = m_data->Stores.at("RecoEvent")->Get("EventFlagApplied",fEventStatusApplied); + auto get_flags = m_data->Stores.at("RecoEvent")->Get("EventFlagged",fEventStatusFlagged); + if(!get_flagsapp || !get_flags) { + Log("PhaseITreeMaker tool: No Event status applied or flagged bitmask!!", v_error, verbosity); + return false; + } + // check if event passes the cut + if((fEventStatusFlagged) != 0) { + Log("PhaseIITreeMaker Tool: Event was flagged with one of the active cuts.",v_debug, verbosity); + return true; + } + + if (useSimpleDir) { + Log("Using simple direction", v_debug, verbosity); + } + + if (useTrueDir && useMRDTrack) { + Log("Unable to use two directions; defaulting to true direction", v_debug, verbosity); //will change to not use true, to be more general, once MRD track usage is tested and acceptable + useMRDTrack = 0; + } + if (multiGrid) { + Log("Using MultiGrid(3)", v_debug, verbosity); //Currently using only three. May make variable later + std::cout << "Using MultiGrid(3)\n"; + } + + this->FindCenter(); + if (multiGrid) { + for (int i = 0; i < 3; i++) { + Center.push_back(vSeedVtxList->at(centerIndex[i]).GetPosition()); + } + } + else { Center.push_back(vSeedVtxList->at(centerIndex[0]).GetPosition()); } + this->GenerateFineGrid(); + for (int i = 0; i < 3; i++) { + std::cout << "Center " << i << ": " << Center.at(i).X() << ", " << Center.at(i).Y() << ", " << Center.at(i).Z() << endl; + } + m_data->Stores.at("RecoEvent")->Set("vSeedVtxList", vSeedVtxList, true); + Center.clear(); + + return true; +} + + +bool VtxSeedFineGrid::Finalise(){ + Log("VtxSeedFineGrid exitting", v_debug, verbosity); + return true; +} + +Position VtxSeedFineGrid::FindCenter() { + double recoVtxX, recoVtxY, recoVtxZ, recoVtxT, recoDirX, recoDirY, recoDirZ; + double trueVtxX, trueVtxY, trueVtxZ, trueVtxT, trueDirX, trueDirY, trueDirZ; + double seedX, seedY, seedZ, seedT, seedDirX, seedDirY, seedDirZ; + double peakX, peakY, peakZ; + double bestFOM[3]; + + double ConeAngle = Parameters::CherenkovAngle(); + + // Get true Vertex information + Position vtxPos = fTrueVertex->GetPosition(); + Direction vtxDir = fTrueVertex->GetDirection(); + trueVtxX = vtxPos.X(); + trueVtxY = vtxPos.Y(); + trueVtxZ = vtxPos.Z(); + trueVtxT = fTrueVertex->GetTime(); + trueDirX = vtxDir.X(); + trueDirY = vtxDir.Y(); + trueDirZ = vtxDir.Z(); + peakX = trueVtxX; + peakY = trueVtxY; + peakZ = trueVtxZ; + bestFOM[0] = 0; bestFOM[1] = 0; bestFOM[2] = 0; + centerIndex[0] = 0; centerIndex[1] = 0; centerIndex[2] = 0; + RecoVertex iSeed; + RecoVertex thisCenterSeed; + + if (verbosity > 0) cout << "True vertex = (" << trueVtxX << ", " << trueVtxY << ", " << trueVtxZ << ", " << trueVtxT << ", " << trueDirX << ", " << trueDirY << ", " << trueDirZ << ")" << endl; + + FoMCalculator myFoMCalculator; + VertexGeometry* myvtxgeo = VertexGeometry::Instance(); + myvtxgeo->LoadDigits(fDigitList); + myFoMCalculator.LoadVertexGeometry(myvtxgeo); //Load vertex geometry + + // fom at true vertex position + double fom = -999.999 * 100; + double timefom = -999.999 * 100; + double conefom = -999.999 * 100; + myvtxgeo->CalcExtendedResiduals(trueVtxX, trueVtxY, trueVtxZ, 0.0, trueDirX, trueDirY, trueDirZ); + myFoMCalculator.TimePropertiesLnL(trueVtxT, timefom); +// myFoMCalculator.ConePropertiesFoM(ConeAngle, conefom); + fom = /*timefom * 0.5 + */conefom; // * 0.5; + if (verbosity > 0) cout << "VtxSeedFineGrid Tool: " << "FOM at true vertex = " << fom << endl; + + for (int m = 0; m < vSeedVtxList->size(); m++) { + RecoVertex* tempVertex = 0; + iSeed = vSeedVtxList->at(m); + seedX = iSeed.GetPosition().X(); + seedY = iSeed.GetPosition().Y(); + seedZ = iSeed.GetPosition().Z(); + seedT = trueVtxT; //Jingbo: should use median T + if (useTrueDir) { + seedDirX = trueDirX; + seedDirY = trueDirY; + seedDirZ = trueDirZ; + } + else if (useMRDTrack) { + iSeed.SetDirection(this->findDirectionMRD()); + seedDirY = iSeed.GetDirection().Y(); + seedDirZ = iSeed.GetDirection().Z(); + } + else if (useSimpleDir) { + tempVertex = this->FindSimpleDirection(&iSeed); + seedDirX = tempVertex->GetDirection().X(); + seedDirY = tempVertex->GetDirection().Y(); + seedDirZ = tempVertex->GetDirection().Z(); + delete tempVertex; tempVertex = 0; + } + /* else if (usePastResolution) { + TFile *f1 = new TFile(InputFile); + TTree *t1 = (TTree*)f1->Get("phaseIITriggerTree"); + t1->Draw("((trueAngle*TMath::Pi()/180)-MRDTrackAngle)>>hs1", "abs(deltaVtxR)<2000"); + TH1D *hres = (TH1D*)gDirectory->Get("hs1"); + double smear = hres->Random(); + iSeed.GetDirection()->SetPhi((smear)+findDirectionMRD().GetPhi()); + iSeed.GetDirection()->SetTheta(findDirectionMRD().GetTheta()); + }*/ + if (useDirectionGrid) { + for (int l = 0; l < 50; l++) { + double theta = (6 * TMath::Pi() / 50) * l; + double phi = (TMath::Pi() / 200) * l; + seedDirX = sin(phi)*cos(theta); + seedDirY = sin(phi)*sin(theta); + seedDirZ = cos(phi); + myvtxgeo->CalcExtendedResiduals(seedX, seedY, seedZ, seedT, seedDirX, seedDirY, seedDirZ); + int nhits = myvtxgeo->GetNDigits(); + double meantime = myFoMCalculator.FindSimpleTimeProperties(ConeAngle); + Double_t fom = -999.999 * 100; + double timefom = -999.999 * 100; + double conefom = -999.999 * 100; + double coneAngle = 42.0; + myFoMCalculator.TimePropertiesLnL(meantime, timefom); + myFoMCalculator.ConePropertiesFoM(coneAngle, conefom); + fom = timefom * 0.5 + conefom * 0.5; + + if (fom > bestFOM[0]) { + if (multiGrid) { + bestFOM[2] = bestFOM[1]; + bestFOM[1] = bestFOM[0]; + centerIndex[2] = centerIndex[1]; + centerIndex[1] = centerIndex[0]; + } + bestFOM[0] = fom; + peakX = seedX; + peakY = seedY; + peakZ = seedZ; + SeedDir = iSeed.GetDirection(); + thisCenterSeed = iSeed; + centerIndex[0] = m; + } + } + } + else { + + myvtxgeo->CalcExtendedResiduals(seedX, seedY, seedZ, seedT, seedDirX, seedDirY, seedDirZ); + int nhits = myvtxgeo->GetNDigits(); + double meantime = myFoMCalculator.FindSimpleTimeProperties(ConeAngle); + Double_t fom = -999.999 * 100; + double timefom = -999.999 * 100; + double conefom = -999.999 * 100; + double coneAngle = 42.0; + myFoMCalculator.TimePropertiesLnL(meantime, timefom); + myFoMCalculator.ConePropertiesFoM(coneAngle, conefom); + fom = timefom * 0.5 + conefom * 0.5; + std::cout << "seed (" << seedX<<","<Get("MRDTracks", Tracks); + m_data->Stores["MRDTracks"]->Get("NumMrdTracks", numtracksinev); + + if (numtracksinev > 1) Log("Multiple tracks need work; just using first for now", v_debug, verbosity); + double gradx, grady, theta, phi; + Direction startVertex, endVertex, result; + BoostStore* thisTrack = &(Tracks->at(0)); + + thisTrack->Get("VTrackGradient", gradx); + thisTrack->Get("HTrackGradient", grady); + theta = atan(grady / gradx); + phi = asin(pow((gradx*gradx + grady * grady), 0.5)); + /*TRandom3 smear; + Direction vtxDir = fTrueVertex->GetDirection(); + Direction result; + result.SetTheta(smear.Gaus(vtxDir.GetTheta(), 0.4)); + result.SetPhi(smear.Gaus(vtxDir.GetPhi(), 0.4));*/ + result.SetTheta(theta); + result.SetPhi(phi); + + return result; +} + +RecoVertex* VtxSeedFineGrid::FindSimpleDirection(RecoVertex* myVertex) { + + /// get vertex position + double vtxX = myVertex->GetPosition().X(); + double vtxY = myVertex->GetPosition().Y(); + double vtxZ = myVertex->GetPosition().Z(); + double vtxTime = myVertex->GetTime(); + + std::cout<<"Simple Direction Input Position: (" << vtxX << "," << vtxY << "," << vtxZ << ")\n"; + // current status + // ============== + int status = myVertex->GetStatus(); + + /// loop over digits + /// ================ + double Swx = 0.0; + double Swy = 0.0; + double Swz = 0.0; + double Sw = 0.0; + double digitq = 0.; + double dx, dy, dz, ds, px, py, pz, q; + + RecoDigit digit; + for (int idigit = 0; idigit < fDigitList->size(); idigit++) { + digit = fDigitList->at(idigit); + if (digit.GetFilterStatus()) { + q = digit.GetCalCharge(); + dx = digit.GetPosition().X() - vtxX; + dy = digit.GetPosition().Y() - vtxY; + dz = digit.GetPosition().Z() - vtxZ; + ds = sqrt(dx*dx + dy * dy + dz * dz); + px = dx / ds; + py = dx / ds; + pz = dz / ds; + Swx += q * px; + Swy += q * py; + Swz += q * pz; + Sw += q; + } + } + + /// average direction + /// ================= + double dirX = 0.0; + double dirY = 0.0; + double dirZ = 0.0; + + int itr = 0; + bool pass = 0; + double fom = 0.0; + + if (Sw > 0.0) { + double qx = Swx / Sw; + double qy = Swy / Sw; + double qz = Swz / Sw; + double qs = sqrt(qx*qx + qy * qy + qz * qz); + + dirX = qx / qs; + dirY = qy / qs; + pass = 1; + } + + // set vertex and direction + // ======================== + RecoVertex* newVertex = new RecoVertex(); // Note: pointer must be deleted by the invoker + + if (pass) { + newVertex->SetVertex(vtxX, vtxY, vtxZ, vtxTime); + newVertex->SetDirection(dirX, dirY, dirZ); + newVertex->SetFOM(fom, itr, pass); + } + + // set status + // ========== + if (!pass) status |= RecoVertex::kFailSimpleDirection; + newVertex->SetStatus(status); + + // return vertex + // ============= + return newVertex; +} + +/*double VtxSeedFineGrid::GetMedianSeedTime(Position pos) { + double digitx, digity, digitz, digittime; + double dx, dy, dz, dr; + double fC, fN; + double seedtime; + int fThisDigit; + std::vector extraptimes; + for (int entry = 0; entry < vSeedDigitList.size(); entry++) { + fThisDigit = vSeedDigitList.at(entry); + digitx = fDigitList->at(fThisDigit).GetPosition().X(); + digity = fDigitList->at(fThisDigit).GetPosition().Y(); + digitz = fDigitList->at(fThisDigit).GetPosition().Z(); + digittime = fDigitList->at(fThisDigit).GetCalTime(); + //Now, find distance to seed position + dx = digitx - pos.X(); + dy = digity - pos.Y(); + dz = digitz - pos.Z(); + dr = sqrt(pow(dx, 2) + pow(dy, 2) + pow(dz, 2)); + + //Back calculate to the vertex time using speed of light in H20 + //Very rough estimate; ignores muon path before Cherenkov production + //TODO: add charge weighting? Kinda like CalcSimpleVertex? + fC = Parameters::SpeedOfLight(); + fN = Parameters::Index0(); + seedtime = digittime - (dr / (fC / fN)); + extraptimes.push_back(seedtime); + } + //return the median of the extrapolated vertex times + size_t median_index = extraptimes.size() / 2; + std::nth_element(extraptimes.begin(), extraptimes.begin() + median_index, extraptimes.end()); + return extraptimes[median_index]; +}*/ diff --git a/UserTools/VtxSeedFineGrid/VtxSeedFineGrid.h b/UserTools/VtxSeedFineGrid/VtxSeedFineGrid.h new file mode 100755 index 000000000..cbe7770a9 --- /dev/null +++ b/UserTools/VtxSeedFineGrid/VtxSeedFineGrid.h @@ -0,0 +1,85 @@ +#ifndef VtxSeedFineGrid_H +#define VtxSeedFineGrid_H + +#include +#include + +#include "Tool.h" +#include "ANNIEGeometry.h" +#include "Parameters.h" +#include "TMath.h" +#include "TRandom.h" +#include "TRandom3.h" +#include "FoMCalculator.h" +#include "VertexGeometry.h" + + +/** + * \class VtxSeedFineGrid + * + * This is a blank template for a Tool used by the script to generate a new custom tool. Please fill out the description and author information. +* +* $Author: F. Lemmons $ +* $Date: 2021/11/16 10:41:00 $ +* Contact: flemmons@fnal.gov +*/ +class VtxSeedFineGrid: public Tool { + + + public: + + VtxSeedFineGrid(); ///< Simple constructor + bool Initialise(std::string configfile,DataModel &data); ///< Initialise Function for setting up Tool resources. @param configfile The path and name of the dynamic configuration file to read in. @param data A reference to the transient data class used to pass information between Tools. + bool Execute(); ///< Execute function used to perform Tool purpose. + bool Finalise(); ///< Finalise function used to clean up resources. + + + private: + Position FindCenter(); + + void GenerateFineGrid(); + + Direction findDirectionMRD(); + RecoVertex* FindSimpleDirection(RecoVertex* myVertex); + +// double GetMedianSeedTime(Position pos); + + + std::vector* vSeedVtxList = nullptr; + std::vector vSeedDigitList; + std::vector* fDigitList = nullptr; + + std::vector* SeedGridList = nullptr; + RecoVertex* fTrueVertex = 0; + std::vector Center; + Direction SeedDir; + int fThisDigit = 0; + int fSeedType = 2; + int centerIndex[3]; + int verbosity = -1; + int v_error = 0; + int v_warning = 1; + int v_message = 2; + int v_debug = 3; + bool useTrueDir = 1; + bool useSimpleDir = 0; + bool useMRDTrack = 0; + bool usePastResolution = 0; + bool useDirectionGrid = 0; + bool multiGrid = 0; + + // \brief Event Status flag masks + int fEventStatusApplied; + int fEventStatusFlagged; + + std::string InputFile =" "; + + std::vector* theMrdTracks; // the MRD tracks + int numtracksinev; + + std::string logmessage; + +}; + + +#endif diff --git a/UserTools/VtxSeedGenerator/README.md b/UserTools/VtxSeedGenerator/README.md old mode 100644 new mode 100755 diff --git a/UserTools/VtxSeedGenerator/VtxSeedGenerator.cpp b/UserTools/VtxSeedGenerator/VtxSeedGenerator.cpp old mode 100644 new mode 100755 index f7f0f8518..f64cdce59 --- a/UserTools/VtxSeedGenerator/VtxSeedGenerator.cpp +++ b/UserTools/VtxSeedGenerator/VtxSeedGenerator.cpp @@ -78,10 +78,10 @@ bool VtxSeedGenerator::Execute(){ } // Generate vertex candidates and push to "RecoEvent" store - if (UseSeedGrid){ - Log("VtxSeedGenerator Tool: Generating seed grid",v_debug,verbosity); - this->GenerateSeedGrid(fNumSeeds); - } else { + if (UseSeedGrid) { + Log("VtxSeedGenerator Tool: Generating seed grid", v_debug, verbosity); + this->GenerateSeedGrid(fNumSeeds); + }else { Log("VtxSeedGenerator Tool: Generating quadfitter seeds",v_debug,verbosity); this->GenerateVertexSeeds(fNumSeeds); } diff --git a/UserTools/VtxSeedGenerator/VtxSeedGenerator.h b/UserTools/VtxSeedGenerator/VtxSeedGenerator.h old mode 100644 new mode 100755 index 2db99f44f..1da791550 --- a/UserTools/VtxSeedGenerator/VtxSeedGenerator.h +++ b/UserTools/VtxSeedGenerator/VtxSeedGenerator.h @@ -9,6 +9,8 @@ #include "Parameters.h" #include "TMath.h" #include "TRandom.h" +#include "FoMCalculator.h" +#include "VertexGeometry.h" class VtxSeedGenerator: public Tool { @@ -65,6 +67,12 @@ class VtxSeedGenerator: public Tool { /// Use VertexGeometry to find the seeds /// \param[in] bool savetodisk: save object to disk if savetodisk=true void PushVertexSeeds(bool savetodisk); + + Position FindCenter(); + + void GenerateFineGrid(int NSeeds, Position Center); + + RecoVertex* FindSimpleDirection(RecoVertex* myVertex); /// Data variables @@ -79,10 +87,12 @@ class VtxSeedGenerator: public Tool { int fLastEntry; int fCounter; int fSeedType; + int fFineGrid; std::vector* vSeedVtxList = nullptr; std::vector vSeedDigitList; ///< a vector thats stores the index of the digits used to calculate the seeds std::vector* fDigitList=nullptr; + // Initialize the list that grid vertices will go to int UseSeedGrid=0; std::vector* SeedGridList = nullptr; diff --git a/VGCheck-pdfcharge.root b/VGCheck-pdfcharge.root new file mode 100644 index 0000000000000000000000000000000000000000..315a9464e2003f1deab8e44ef2f72878e0eff577 GIT binary patch literal 19920 zcmeHv1ymf}nr-7Qfe_r15Q4kAOVHpB!Gl|HcY+0XcMSx$27(1AxNCp_p>dkVUc-MU zcVx|-d+*Gw_11eo)m3#)HK$J1K6`)P(RQ@6a|QrT`v3rd2>>8J4&A+>VHM~O1Kl0r zp$Ag{KzI@WKraKpj8ICRMZbM=7#s= zX=!ar%_MH*Y{aBuX=`q6>TK%d%%m2 ztpEUABQzsN=>C`^03bl}U%dv5_rudthCilPcz9YEoQL%As?sdtPppmX?M+PGoK0;_ zOii2}O`S*{uTV(*qi;=aq9P%p|KsY794d&N`Kt#9FQ-n3+Je==g za|%zo%rdCJmqc7FbMPk-ewcZL19BWAQyWA}VNEEFIqH*gbSwDzl99Wxqi9=tc5$A= z_`wh+YlMA_Xjtvn1C&FE8GGa8F>s_(kQd`3>=Tp81)JxSN%+Ape#r=-3*P4PA9mKR zj3Ac|eKr|Zr`4eSdG13ov?5~;uQ6ILM<8269;zbU>i=A6jv7C`vUlN^Dq{7o`{`kH z=6YSr!JLs=f0NQf2cWC^mA?TN+(_q}oSvfcOV_d!*FUMl0$ zn?xluQj%8F*rbFV9L@|0>bj-QjC5E#YE@6Qy)cVM?NMtwLeMna1;r)x%)z>|oCB2^ zBxxjsu9Qo^)VXJmD-)NduWyWZ)~mWQLp$1hGB^fCDZkLML#7TsR1T3p=@Y#D4hMRQ zIZSj%K+v~0-Pmbq2`NNcVITpY;xhUgHelQzO0%_QN@RFnz28}|A-=udA%`Fu9M}qy z4qIb3h+l`g+D7Tnd>y)CYvifD0+IaJq?xldZIpj_2(w z2hO3)_KGN|g=iWRDz=Xy`Z z=R0-o3%Q@c`SRqvv#utQ$$@8)IWCP7xk*)>Zu?|FtRTNkA+`%VTXVao5 zHG(~*Qeu0jx9N8uz2FnvNwxQ*7*?r$wzc1yZwh`VH^`dx#w{jy--vCS27eRVh;@r> z098XuLM5E$nVe2|Pra$^!n3)Fa2j74i0|QyV5#hv4Kq`u&-6FwB zZx*T5L&vt=eTW+utU52kEjMcJwv4RFoWAJ*r)TW;-?@89)SbTb@fwG`bJ!Q_zp3GN zk+_+JHt2SY02*%Kso%ME>7m zMYgxIw0$s;xJnO!|6mvG-~g_F$n%YzJXD^If?Gw1wBpj}j4Z+O7Ym7Ari+1UM?X@hCW+?djhl}QRy*<4L48B*(wGZ9EffAsK^dTqcU$9!kExw{T#cI>xJ?baE zgLGOGpA-Qz1x6J!z_ljV1B3V{7vP`$*Dd;W*;!-OFxy>ATqs6byEE&Vy_q}QwPWX- z8BM(E#3h_cHaQEOzMo<4F1n@?zkMW(XwpC^<;k0sH!#&T7LW0Td`Zr^)N1Q(nsn=I zntfe)LlK(uvp?r&XjJh_MPzv}#*ESy??2;8WMJt~SHPkLX?Xz^ zjUlV)6K-W1JoD_BwEbpLJOS+2?6#P6)gR+-+2w$|Ag!6^)2UIs89px<-H+RisUBGDnD{A0 zMjCg?;&XycnpeDgeODFZp~;zi%UWFGw5oBE1)d0}Y@xnleM?P3QYxiz)uuptNnM1+ zBCa@b?IlO*5dro|>`*Mj%pP;{t-jWjJMw^`!WQ9r|Mzgdu5haW zTZR|(P|Z3#YS!#gvvd!dMLGJdSp*YPYiFZBj?k)$035*VFCO%bnlg5O49Jv!$wcTG zB1464fdqjK?3)F89s(JFb*ZA=cqX! zre7^tpV|LHs^{Ty%RbG#8XazT$}VyZ?bn@dl^vOGbF!(Txk1yIZq`1HI&PiJKOuDF z&?rU!ab>^qi?-{OCquiNrU$aes5^a*oodgp?sQ+IQ04x?`Qjo0Z5^|5sJX~4%0y#q zus_MkGlXx<$NiSwE;d`N2Y4QXUvN0Xg_Gz8BufZSOEZ#L3)3#+^PXM+{Z)UuT` zQIZ$fR~U4?^1WW^YN|1!?Ype&mfH*Bu3UDfa;%@}tZ&tRN26w4ujE*lXChIHHqEQ1 zk=o0ve$8Zc6SvnTxo`zDeQ+J{zQmUw0&{+8byxTPWVuX{#(hlSl>+H1QU?+^rj)#) zJf!x+4wfZ}w{X3Wkay?=*z|T(rQhpjxCRn-VT|d17{60oJ(t-(rn>>AaGKV*@jB)| z>q7cAxf*u|^+Mjc;ZEeS*{7nL8he{O)kt0vOW;Z@AqSt;oNr>pF_s%ThuJ$H9M|)O zrs~9?dh=}X{0xe}>?EsWIX`1a53Np7G+=m1zAu*H!CufY2mA@_l@ivP%B^3Si}2ns zucC6{JbZC`jOV`0Y`4WEV1;z*}$#D zV9WqQENI;oruI>AE{}raeGnXm(Qm;a*_ayHLVYBTJoNk@T_g;^xd*jn_qcltN z2V7w<`F57lua|FXb&IZ4YEvt;98)a9WlyB%#78xsg8R)y%v)&ovO<;DWbw%K!BVwT zltb@}hp5P@-srn0RXk5UBNB-FeEmJVR*Qb@fRvUmlza}8| zzoR>rUtDZIj$hQ%&(>*;Chy+rW7~J+kZ>@}qoFn5zUb$|S1-7)Hz@fO06_M8{M@j> zUg-rd9~ickw3M~vFYKiTUE->}=*gNvvL;$$6>TeaQL$B6r{H+gI5q<@R^Xjbvw8e4ySYl3(CC2?wVpzDpOAOK4+0n?@>dE7g!rsW)g2LL= z*4)_wYAS%z6{$B(gK_2S~E019WziU+vl%LGT18f|jA`Wn|0*EAYfZC>Lm>Rp{X(2>XX z7ro^Q?0E9bE>-Vdvxs-MbElOoW$eZ&hv%8qNwH+$w?Hhfntgv?;QBpzN^4g(_4gZI z3586z`)|Hi>-|&ZnWu$JI{QEOd)uNC)X#)cEP_I|j$}cZ*Ssen7ab0uySbQtbTBD% zZJe9LH_kPM=aD&ydRrPcsoW7Zopf3D80feYde7eba^WtEbnJ&mzh#EeD!DvhiTB#g z9~Etvq7GCDqC9&~&np}3u^LV7!SUq1Z$kJ?^PY;u8G6 zMtn~iIzPPPH8i4o##q>GsyG*Lfjfbv*#ysMo_DKBpozU2IODUL+%nePiMtXL%&WL# z;lq;Xv4)~HRGm22Xt}W>kegnL(FC*CVMc#dgm4VI=bF3Ig}CUq(-qv5Q~5GWtV6kE zOmx*~P{VGtBC)5ji>Co*)vp6#v?mUTS)16ym0OL+9SC7RnHCtRR^;eu@K=V z5(K1yk&#`+vt$F>f<-(-2vaD(($xxUoFk#LOx(V-m0JSi^4`JW_c6SoVq79{?c4Jl zBlX92*eQxA0P^mvl$g&lVD&gYl6(HX4NTrqh)O4wAYpH2*l>MYNWS77; z{e5sla=wvMnU*8I?)pku_KM2jwqky)ZT8nd%vrQ+Qt7;bB4{Cni0d>xv3TnjnsZRe zdh$~pz7;Rq4kv!s=F|SrN|Qj7L%K?ze5Rl-9iM4~kZLmLwJ%~=au-nA3q8xfZXLmV zMeenBDu%Um9dz$w7>DS45D~vPKW5!O@1LEGxN5rE0Q5Ni<-*{T`etF^^}Ea9aD}(| zWLZ^0r&ZBguFEVoS_D@?>U-Ayg20&MoG6>|*{<|8j?biRsM_~<7W+ViOazKdmF*bECQ86KOq7n!Cd>cs}3^C<&XhB!CEJJ^q1`=G>+i{mx5gM<8MZPFXH zmeh#JR)cQ~>aNSfUN-E6vuZc)&{~&A!{52*epEnjxIOu_m*MAzS8;z>an4dq654R- z1>fNs1~m37@{7inC;kUT|Vxx&56K*TIC-$-M31Ht z7)P!GL)QFc3Kdm85)nB7)Bzwi&J9m=a5cMdY=b+za0GGPr_QI&dt>GLV0I(t>W%vW zFsF92^^N>OT559Fw+4PlhTxJb=n{VC!N5kbNl2Za(iGA&yqFCg3)j#yFR#cCA~|=6I3l;03zPRs<=0@JhS6*qcw6 z1RVq9ZDsP#`~j`G>ISv?>l?^agpl^`$Z78LAUQABq^lnYP7@NwYuy>otDK*+T+dK3 zGI{oZ0wZ{A08RYiJ;QI{IuD#eQ$B)wZr(*|YAW^UE>UgUu{i}3{5+7)PSJAJA{X7R z=)W1EKe&`Je~-t)cDwlnlBEQD zcZ1Q_2Ke!t^78Dc)qbzQDBmCzZMRuervTS-+BT zlxbC_{yM<@9)`2`P!US^M_+(6{z!L@2f7pdwf`H*$ivdc#n}SmFX!KcUxfjTKM?-c z+jfNnk*5L&5!ReU`5k_goa4&ZUxNdXG`|R~9`-Zp5KuiK_6<-~SsEc1-v=QI(zR#Eg_>GCR8t|3NGV zxRk5-kzOfJV4C7#0$R+Nua7LsckN6Gz%^t{Yz4_c9x2slEtC#R?k6qM_BFN5m229) zpb5mEPlL?uXUM%*4ktuaA*)X{?uHvnO;<}(2jT&~of%A?-kknwXLan9YK_fSssOvx zTfwO%<^5KaN_17tfW>Dxt}bBxi^bfte)QN6?E2gwq)oB@to<{-iJPSTWy1`QF?S6) zYP(iJ`*z|HkhveT?Nd4Oa3emr0#F-8h_&2BjqNtE8F4p1TnxGRq6YU)hfLe`J`Z}g zXF?nTVau3DYK+IaAV|hU&-61tacf8*x(CoA!^MVsM{wpVZ*$CfjgS`>L|_|#Xv{;h z%Z(7@Bf7_`VH{npMSj zl9;LJF}AaXwqu|c5cMJSn;T~Uz&$*a8%~CX zS0kMVot-ZSI@dcpa|Sw*gM>SWJ3G5lGqRJ@Qj)Xc2NB-U_@xR?a%v2Peh|hadY?k~Pl1{&~&wFDWpuE4n+oJ3GNx^rfclVBTZ~ z-^ikoS+QFRIAQCY$Y;>o zjEl09}~G`vQ9d8pZ{wgHG0WPUOnJW0Xxbgi?KVxSWh-Gst0 zgbzjrg;}4F!i)kXnza-c_VaFjv}Jx#EHNzLC8tzDTe?*`jLHwYNrBu6juRpb%1ud<5t#Vxs zkH~bh-zJ{OHGfex%sXX&qHpvg9F>|#hRC2`ceT7MDg90W?llsxm*I}xKE z>laYxL>#m=n3qd*DXutA#8B^peTVd`4Z>0E{KT1@#ILj#Gm4yMHb>{vc3+XovO*$X z#lcyk&FA$TgdUW+^F3 zs_-t)5!fN!nsFueAk>vDG;AG%XsoS@uA;i(())9DB21P(ZSxA6X+-vgVb7L%HWv6M zXGBM?mtZ?NcJVejH_B<;j=B(xK|t^|2g*OV2mT@3{Z+7}Yrpx2^2oV|5lhrR*!K{G za<2~tV2Au*%oACxRG@wH_cBE8W|pF)B)o0!UUS4Vs!P{X&nQo=QKX<_6igR~ykt+Q zs{?<}D~DjfGDiCZg{q(+JCt#oS~oHVFm5zoX5i~v6pZ{~+NXz+Q$`m<9+N0}%_d!;f8xsC1SL&F7OI11b{ z!%4zTioYLoc7AxKB76|DJp4n5tiV=2xtj6z=46z_&~fD(%p6v2n`7+TDIi_uvi<2& z$aaaen5Xw919|%Sy@)y25^v*rDRK02ERwhMsx8uXyv5PxOJ5-UuzJXA-7Cnx%})_7 zuo9;!vZaJGB}L#%tkf4QkRWM<`>I{3vzM82hOGYI;QP`5mUwz`K~;1b39bp3-s$5n zaFP=D{VnWXTb~7us2TAwED@qM!sV{weh@-&McNB($@P(Hn#KJ9EXTY&Z|9nw0+u=l z1Z?>K4yOjhzCc)fjJJs|6PLx>o)yvlwWKPe%m z5*|^H?Vbz?NnOP+KodNUK=yJk(O#bwY6lkx_RyzT( zmYC@1&FH$qZwbCxBhEzhL=3+}%=wH8hXx0?DfGD@>3m~fx8a5*jne?@JP;jFgDdo}=dJlH0N2#BoM%wIDB8zL}Jll|Q=L)n8JUQ*1@9vJ{&s}aJ3A|>;P zlfF0fue*215o`JlnftG!ZCh8Ny+#(N;Cruf5m7Wpyv39 zkt*CHMx_4)3EL;Iw+~|<7@6Ifl z<~*2-I(Yb0bWFU=EQke>Fp_~h13_?bffR6b%(4D{gop4i@jBg1rNuU4SVNVOxid82 zuKRN50HXM$e%}ovx^kNU-Aq}4TsDf0{Fisw0Kh9W1yyLJpbrWpZd>#tkRlI2K3D$> zkgVcQv_3eAKR7;B3sYw!10!o|EG>!;P85H=1_jiM2oT+d3pGItieVY%q)11N>EcNg zb%mbDM9w9JJp`g@Zxh}x+fhW!k#TkjS>u#i5pR6Qs4JT4d5-BqC4s&8CD~6OzN2$S13pl8>KLNgiu>x$l{awA9PV4lR<6)(+g$p) z!`9x8yZ!iHuZRr`2ZMW!dn-ldVWFm*!NT}*e077nzRNLS+sVP;!ouxPrUs*nVB72} zi^1?nCgH)Jfp101KIlgA{=M&gnXg@q&ePr{EX51ZP1D?}@9l8XRd2@-U9tv~z|nfs z=={!@0a(Fz{zdOfUK6N&q&E}@-uFH5-Me$jzWvVUd;0^bua^5Qd*5E-_UU)0bd4)^ zK9Kv_r@e_Iu-P5?orxm;F?P>(>0L(!T^%gIh6gF@_*`GfTVC8q1s>(CTK7y|aIHeL zmzum?sxJlXb}E@J5fwo`pxe_+pxe5yxsdw?HL(5Sk_#X6e1*S}XY1m6JG7FCCIacB zH`muoFD0X)h~`C?lht-#uLA>L&&!_b703ugL+}^`^m_lp$#(^K(b@&_W7yZjp%C^1 zz`(*IWB(cW6NLXgCJ0Bt@@N`2S^Nvr{0HH0+N)XLi!35Rf` zq}KW@N`q=H;NI7rnljuU)0DASyl+kam5m*T_jgVk7e)G$h<|#W|DMMULD@q4n&qM1 zAgtbE8e=`x*5*<1#QS7{oM^&$;my`G@uW|WWCl}lOVn;Iygw`a+dgbIWYw3kNVKRp zVt;@?+$KzxEF#vpt#v3!oUf7^%l1!2oaZ7G_k?7 zlVxLr{(EdxI27Xw{<=2mpqn$tFQ}H2wEpay*-YA|MfCn!PmY|1W{W2M?#2~HksM2N2S!a{m}O# z)wKK?o@x>{>G9d7i`TYEFJkR31_`I)Xk6}j5~Q?Y>xa^`=6(Pcf7?vjCj;-PSVaFc zf`3-b{}(IfYhr{PRonav$y#igyMSTt1c{}qYZm|_FZvA3iJn$-%;Xw)PyM!tO+*=4CSLBTqGTO35Q88eXR zaoTUUQkCgmvsl&5XPwU0PI~y-)jjKH65r?bNyoO8fS*(825OKENi~_povBQqH=S!b zV{96lS;sc)2Gp{SEBWgxh^3cG{Mpj@vjY8pVFkKQkAk|Tkc1K|nYAxA273UPye5Y# z#wjl&PirrUZMTZX-Hg=<6meYh16FH4$~ho|&}LcsYoBHUzHMa9P*NtUNptHYLzH%Q zv8y`u8$KQ2;kofSq5l2fxaOM*1c-7Z-E!L zNsK}Y1fFZc;#q2vHO3qCuW6VU|5GS`3g!Poq5NHi0v2BL##9@jfX<3-d>9Udvwa-J zdR(^*T|Eq43H%SPTZSgV0s#IemlZ=F>VIgREc8zQ^ddfJ+P`-(5%l%{$z>4G84Ul_ zbbRQ;{gZQ=p*jAe6F#9y{_#n3P)z>e$yCsc{;Byy&~*RsJP_#f{qqwCpo#yvvGo7; z*zw~p6+mVC`q`swg&zi=AAY~!LAG+1PR{5m&_$_6Hl~g;wq|xX;&v|YtlcTBEp4q# zO`svapSBXNXq@$ z8P6N)woI^me4 ze^iQ!(gK~f58cdKtw1g0II}HI=4PO!RRQwL45zVrk$g`YMIO-!aDZdENVBV;FwLlx zzo~lVTATTGP@~dfptGqPp0-XBnn6IILch+tUAlsGNS*L%ddTRa+>U)@5g*>rN0P-F*ym}Y0znH6 zXKXKE46*epRSwY?5Wd@czqLPOSIE|Lem+b4Iwy?01PvYglV^p3Ghrp=#zHnamIDX3 zbOaB5_2T%gr$uObVF64^>dE^Mt?uiMzKF0vDhyKV+)Z&=W|`SM1=6M8 zoRVA`bIh6EFTO2t&;7Ecg|~eCJg`KY+J>-Ff~aZ??a+dV%TAF`xe2%>&fBx8e#Xq2 zuU%4oB@!HVk@Ep_pa-aeZ}uH4I28$(p4rLun1xNZV=d^(?V;6(|u#rmh0HXxhWizj`dh`7mL~$XF|zgMO$FMm2swI z)3b52)Xd#VUtwimWzu}&TI>h#Fx7REiMRlH*2f@5<|-=IQ3@9AZJJ>zOZ2hB zIJ@T3%U!ng#uhEq{7!l)4;1w`Os!@FMLR1^zT4zL?{bGRo4_KnFzU-&k2WPB9s6=p z5e>$h_t+(GYzecE-nB2#7m$%HQ)r8ZrdtRhyx5yEHa3ex>>SD)kivP6%==36)DLN> z^tyqJp|~CZRLuSGp2xfXBw3m*Ob*8LH5w%_NqBXg+{vD%$ia^tBlITYSh}@5_u@%HmlH?eBg1Ri zNAleC+zuu8W5Kf`)niwisAIgcpU0ZJQ-7+xpt!P%Y9qMdhp%EmDP*`c>L81{gJUWA zRnqL{I43W{RB+C(OF-$QL??pFuC?PE?5NF@v3o@XKjzCIzvfm=GJKM;JTAysW`8UZ zBwNZisZW*9QdjL%;rg;Q-iT$;Z^Sr)9?@{#R0K>S;XIm-$s!KRB6)cO!rmk-#$o6u zu<}~uAsxPL?UdKYKw7l-opI2EtT7p{t{Y+GUG5!2Am(hBsjdXE^u4#{GG8azYM+q} zIA-@hJp|P_6sIHGZ}ti=hvkS-JP|=HdBejTb6DEX-=gPM$3^)HiIho@lSSLn&%`$OF-N` z&Wh5~FZtyqq~U^tiv7;vYcXs2UrJK*8sxsU29f2%jZ$q}VlCn>M}GeLdOB|=L*?0l zm^4vr)wJ=uA<=_8Y`&yd>alL~kw;YeEIMCHsTOfC4`tnSaGDVOhkDRygt#jXATt(& zS`Ff!T7*Gi?12Ng?DVuZxt#+lK<^ssW?MIBHtW|s`rN!MPFpFAz}@_KHmS3G#mtme zTLgO&tBUj+OQ;3uei9_|)a-+T)3G`ph$%kAs_pAY(4F^1K=WbHS2SKOk6zRv|1UV< zp>9E8caaMXTOU=x0~EBgmSc++jec5BQ+gGNw&OZ=A=1oDG?XC){zYy)f^p9_`~++e zof=5BC%p4sfs>M|oeCXQ>w#i^LLO&7&d1_E}dv zbhKdIm95RK+tb+DA($MCo4&+Ph+=W3=T7WDPn)yxe1+++(-2K0TJG z{dgpuJeQBN?>*VGdOb<3YP8GA)GK|Kovo7SKjF>BJRkLQq3TnhAY}@PT#P|)X3#+T z_0dNuv$tajOAe#i5FyVJ zc1LU3RE6^Mi{*EugZVUSp5mjf&6qMKXZx%#X2Wo;$6o~BzYr=(*O#acR`2wm$5)>B zsAuqwiy>x<{h(Ct@rj*(WJsL#&T~AgYQ6#QsG%Zxw^R?8-FI&(X^M$c$c0U4G z=RK`mY)G-)fNGYsRd4Lbsq{ueJY5EU9D!V;o>TZa9B1??mk@)MWrE> zY4Jd|3QYbPPgd6We$tA9tIFh>-O|lN)Fk~|{@FY}j}rMpqzCvD?>H4+8WT-W!?D7O zehKf2kN8cCnpPZq(Jy#rY5*)IQyv5UcQcFNjLIHZDDlHKY+TL8zA9}){1;HIEa1=y ze`FekCrcZ11UsM&mwvJ=;Sm#&^^)t?1_I(*j?wTWGUa@iY5oieCAzQp590F+zrp1q2$g^@iXiwT=|R~iLe}a1;yHih+e!W%^&0+;Pxh!rM+3Xbl+wpIcn+LWnD1ZZAQZU^3WNx2;rO%nz6cxOpn3s58 zLdyvTJTV zAnTXC&oM>j#Zx+0fVn~>qM+l5-~scrGg zd*GDISRd8?2KV9go+z#^7UkkMxVi~E1piAfZhl}S1BQLA$^EiZ`L z(mI_kJe_fW6}uo#4Mv*