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 000000000..315a9464e Binary files /dev/null and b/VGCheck-pdfcharge.root differ