diff --git a/CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h b/CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h index 0654dd742e63d..6bc776a34d580 100644 --- a/CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h +++ b/CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h @@ -21,11 +21,11 @@ class HcalTimeSlew { class HcalTimeSlewM2Parameters{ public: //M2 Parameters - double tzero; - double slope; - double tmax; + float tzero; + float slope; + float tmax; - HcalTimeSlewM2Parameters(double t0, double m, double tmaximum):tzero(t0), slope(m), tmax(tmaximum){} + HcalTimeSlewM2Parameters(float t0, float m, float tmaximum):tzero(t0), slope(m), tmax(tmaximum){} }; class HcalTimeSlewM3Parameters{ @@ -45,7 +45,7 @@ class HcalTimeSlew { HcalTimeSlew() {}; ~HcalTimeSlew() {} - void addM2ParameterSet(double tzero, double slope, double tmax); + void addM2ParameterSet(float tzero, float slope, float tmax); void addM3ParameterSet(double cap, double tspar0, double tspar1, double tspar2, double tspar0_siPM, double tspar1_siPM, double tspar2_siPM); enum ParaSource { TestStand=0, Data=1, MC=2, HBHE=3 }; @@ -53,7 +53,7 @@ class HcalTimeSlew { /** \brief Returns the amount (ns) by which a pulse of the given number of fC will be delayed by the timeslew effect, for the specified bias setting. */ - double delay(double fC, BiasSetting bias=Medium) const; + float delay(float fC, BiasSetting bias=Medium) const; double delay(double fC, ParaSource source=HBHE, BiasSetting bias=Medium, bool isHPD=true) const; private: diff --git a/CalibCalorimetry/HcalAlgos/src/HcalTimeSlew.cc b/CalibCalorimetry/HcalAlgos/src/HcalTimeSlew.cc index 34e2da8eda6df..5f6f3f3f90f58 100644 --- a/CalibCalorimetry/HcalAlgos/src/HcalTimeSlew.cc +++ b/CalibCalorimetry/HcalAlgos/src/HcalTimeSlew.cc @@ -2,7 +2,7 @@ #include #include -void HcalTimeSlew::addM2ParameterSet(double tzero, double slope, double tmax){ +void HcalTimeSlew::addM2ParameterSet(float tzero, float slope, float tmax){ parametersM2_.emplace_back(tzero,slope,tmax); } @@ -11,8 +11,8 @@ void HcalTimeSlew::addM3ParameterSet(double cap, double tspar0, double tspar1, d } // Used by M2/Simulation -double HcalTimeSlew::delay(double fC, BiasSetting bias) const { - double rawDelay = parametersM2_[bias].tzero + parametersM2_[bias].slope*log(fC); +float HcalTimeSlew::delay(float fC, BiasSetting bias) const { + float rawDelay = parametersM2_[bias].tzero + parametersM2_[bias].slope*std::log(fC); return (rawDelay < 0)?(0):((rawDelay > parametersM2_[bias].tmax)?(parametersM2_[bias].tmax):(rawDelay)); } @@ -24,7 +24,7 @@ double HcalTimeSlew::delay(double fC, ParaSource source, BiasSetting bias, bool } else if(isHPD){ rawDelay = std::fmin( parametersM3_[source].cap, - parametersM3_[source].tspar0+parametersM3_[source].tspar1*log(fC+parametersM3_[source].tspar2 )); + parametersM3_[source].tspar0+parametersM3_[source].tspar1*std::log(fC+parametersM3_[source].tspar2 )); } else{ rawDelay = parametersM3_[source].cap+parametersM3_[source].tspar0_siPM; diff --git a/CalibCalorimetry/HcalPlugins/src/HcalTimeSlewEP.cc b/CalibCalorimetry/HcalPlugins/src/HcalTimeSlewEP.cc index 610f6e6529377..1dbfa821e25b3 100644 --- a/CalibCalorimetry/HcalPlugins/src/HcalTimeSlewEP.cc +++ b/CalibCalorimetry/HcalPlugins/src/HcalTimeSlewEP.cc @@ -60,9 +60,9 @@ HcalTimeSlewEP::produce(const HcalTimeSlewRecord& iRecord){ //loop over the VPSets for(const auto& p_timeslew : p_TimeSlewM2){ - double t0 = p_timeslew.getParameter("tzero"); - double m = p_timeslew.getParameter("slope"); - double tmaximum = p_timeslew.getParameter("tmax"); + float t0 = p_timeslew.getParameter("tzero"); + float m = p_timeslew.getParameter("slope"); + float tmaximum = p_timeslew.getParameter("tmax"); hcalTimeSlew->addM2ParameterSet(t0, m, tmaximum); } diff --git a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h index a5fdcdf61b8b9..a7c583e3ee17f 100644 --- a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h +++ b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h @@ -30,9 +30,6 @@ struct MahiNnlsWorkspace { //holds data samples SampleVector amplitudes; - //holds inverse covariance matrix - SampleMatrix invCovMat; - //holds diagonal noise terms SampleVector noiseTerms; @@ -43,36 +40,19 @@ struct MahiNnlsWorkspace { //varied in time std::array pulseCovArray; - //holds full pulse shape template - std::array pulseShapeArray; - - //holds full pulse shape derivatives - std::array pulseDerivArray; - - //holders for calculating pulse shape & covariance matrices - std::array pulseN; - std::array pulseM; - std::array pulseP; - //holds matrix of pulse shape templates for each BX SamplePulseMatrix pulseMat; //holds matrix of pulse shape derivatives for each BX SamplePulseMatrix pulseDerivMat; - //holds residual vector - PulseVector residuals; - //for FNNLS algorithm unsigned int nP; PulseVector ampVec; - PulseVector ampvecpermtest; - SamplePulseMatrix invcovp; PulseMatrix aTaMat; // A-transpose A (matrix) PulseVector aTbVec; // A-transpose b (vector) - PulseVector updateWork; // w (vector) SampleDecompLLT covDecomp; PulseDecompLDLT pulseDecomp; @@ -154,7 +134,7 @@ class MahiFit FullSampleVector &pulseDeriv, FullSampleMatrix &pulseCov) const; - double calculateArrivalTime() const; + float calculateArrivalTime() const; double calculateChiSq() const; void nnls() const; void resetWorkspace() const; @@ -183,7 +163,7 @@ class MahiFit bool applyTimeSlew_; HcalTimeSlew::BiasSetting slewFlavor_; - double tsDelay1GeV_=0; + float tsDelay1GeV_=0.f; float meanTime_; float timeSigmaHPD_; diff --git a/RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFunctor.h b/RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFunctor.h index b715ab52889b2..a3eff964969f9 100644 --- a/RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFunctor.h +++ b/RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFunctor.h @@ -23,7 +23,8 @@ namespace FitterFuncs{ unsigned int nSamplesToFit); ~PulseShapeFunctor(); - double EvalPulse(const double *pars, const unsigned nPar); + void EvalPulse(const double *pars); + double EvalPulseM2(const double *pars, const unsigned nPar); void setDefaultcntNANinfit(){ cntNANinfit =0; } int getcntNANinfit(){ return cntNANinfit; } @@ -37,6 +38,7 @@ namespace FitterFuncs{ void setinvertpedSig2(double x) { invertpedSig2_ = x; } void setinverttimeSig2(double x) { inverttimeSig2_ = x; } + void singlePulseShapeFuncMahi( const double *x ); double singlePulseShapeFunc( const double *x ); double doublePulseShapeFunc( const double *x ); double triplePulseShapeFunc( const double *x ); @@ -52,6 +54,7 @@ namespace FitterFuncs{ std::vector acc25nsVec, diff25nsItvlVec; std::vector accVarLenIdxZEROVec, diffVarItvlIdxZEROVec; std::vector accVarLenIdxMinusOneVec, diffVarItvlIdxMinusOneVec; + void funcShape(std::array & ntmpbin, const double pulseTime, const double pulseHeight,const double slew); double psFit_x[HcalConst::maxSamples], psFit_y[HcalConst::maxSamples], psFit_erry[HcalConst::maxSamples], psFit_erry2[HcalConst::maxSamples], psFit_slew[HcalConst::maxSamples]; diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index 3ec189c6371b9..e977e38afbc43 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -155,36 +155,35 @@ void MahiFit::doFit(std::array &correctedOutput, int nbx) const { if (dynamicPed_) nnlsWork_.bxs[nnlsWork_.nPulseTot-1] = pedestalBX_; nnlsWork_.pulseMat.setZero(nnlsWork_.tsSize,nnlsWork_.nPulseTot); - nnlsWork_.invcovp.setZero(nnlsWork_.tsSize,nnlsWork_.nPulseTot); - nnlsWork_.ampVec.setZero(nnlsWork_.nPulseTot); - nnlsWork_.ampvecpermtest.setZero(nnlsWork_.nPulseTot); + nnlsWork_.pulseDerivMat.setZero(nnlsWork_.tsSize,nnlsWork_.nPulseTot); + + FullSampleVector pulseShapeArray; + FullSampleVector pulseDerivArray; int offset=0; for (unsigned int iBX=0; iBX &correctedOutput, int nbx) const { if (foundintime) { correctedOutput.at(0) = nnlsWork_.ampVec.coeff(ipulseintime); //charge if (correctedOutput.at(0)!=0) { - double arrivalTime = calculateArrivalTime(); + float arrivalTime = calculateArrivalTime(); correctedOutput.at(1) = arrivalTime; //time } else correctedOutput.at(1) = -9999;//time @@ -213,6 +212,12 @@ void MahiFit::doFit(std::array &correctedOutput, int nbx) const { double MahiFit::minimize() const { + nnlsWork_.invcovp.setZero(nnlsWork_.tsSize,nnlsWork_.nPulseTot); + nnlsWork_.ampVec.setZero(nnlsWork_.nPulseTot); + nnlsWork_.aTaMat.setZero(nnlsWork_.nPulseTot, nnlsWork_.nPulseTot); + nnlsWork_.aTbVec.setZero(nnlsWork_.nPulseTot); + + double oldChiSq=9999; double chiSq=oldChiSq; @@ -251,44 +256,50 @@ void MahiFit::updatePulseShape(double itQ, FullSampleVector &pulseShape, FullSam if(applyTimeSlew_) { if(itQ<=1.0) t0+=tsDelay1GeV_; - else t0+=hcalTimeSlewDelay_->delay(itQ,slewFlavor_); + else t0+=hcalTimeSlewDelay_->delay(float(itQ),slewFlavor_); } - nnlsWork_.pulseN.fill(0); - nnlsWork_.pulseM.fill(0); - nnlsWork_.pulseP.fill(0); + std::array pulseN; + std::array pulseM; + std::array pulseP; + + pulseN.fill(0); + pulseM.fill(0); + pulseP.fill(0); const double xx[4]={t0, 1.0, 0.0, 3}; const double xxm[4]={-nnlsWork_.dt+t0, 1.0, 0.0, 3}; const double xxp[4]={ nnlsWork_.dt+t0, 1.0, 0.0, 3}; - (*pfunctor_)(&xx[0]); - psfPtr_->getPulseShape(nnlsWork_.pulseN); + psfPtr_->singlePulseShapeFuncMahi(&xx[0]); + psfPtr_->getPulseShape(pulseN); - (*pfunctor_)(&xxm[0]); - psfPtr_->getPulseShape(nnlsWork_.pulseM); + psfPtr_->singlePulseShapeFuncMahi(&xxm[0]); + psfPtr_->getPulseShape(pulseM); - (*pfunctor_)(&xxp[0]); - psfPtr_->getPulseShape(nnlsWork_.pulseP); + psfPtr_->singlePulseShapeFuncMahi(&xxp[0]); + psfPtr_->getPulseShape(pulseP); //in the 2018+ case where the sample of interest (SOI) is in TS3, add an extra offset to align //with previous SOI=TS4 case assumed by psfPtr_->getPulseShape() int delta = 4 - nnlsWork_. tsOffset; + auto invDt = 0.25 / nnlsWork_.dt; + for (unsigned int iTS=0; iTStimeLimit_) ? timeLimit_ : ((t<-timeLimit_) ? -timeLimit_ : t); @@ -352,6 +353,13 @@ double MahiFit::calculateArrivalTime() const { void MahiFit::nnls() const { const unsigned int npulse = nnlsWork_.nPulseTot; + const unsigned int nsamples = nnlsWork_.tsSize; + + PulseVector updateWork; + updateWork.setZero(npulse); + + PulseVector ampvecpermtest; + ampvecpermtest.setZero(npulse); nnlsWork_.invcovp = nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat); nnlsWork_.aTaMat = nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.invcovp); @@ -366,14 +374,14 @@ void MahiFit::nnls() const { while (true) { if (iter>0 || nnlsWork_.nP==0) { - if ( nnlsWork_.nP==std::min(npulse, nnlsWork_.tsSize)) break; + if ( nnlsWork_.nP==std::min(npulse, nsamples)) break; const unsigned int nActive = npulse - nnlsWork_.nP; - nnlsWork_.updateWork = nnlsWork_.aTbVec - nnlsWork_.aTaMat*nnlsWork_.ampVec; + updateWork = nnlsWork_.aTbVec - nnlsWork_.aTaMat*nnlsWork_.ampVec; Index idxwmaxprev = idxwmax; double wmaxprev = wmax; - wmax = nnlsWork_.updateWork.tail(nActive).maxCoeff(&idxwmax); + wmax = updateWork.tail(nActive).maxCoeff(&idxwmax); if (wmax 0); + positive &= (ampvecpermtest(i) > 0); if (positive) { - nnlsWork_.ampVec.head(nnlsWork_.nP) = nnlsWork_.ampvecpermtest.head(nnlsWork_.nP); + nnlsWork_.ampVec.head(nnlsWork_.nP) = ampvecpermtest.head(nnlsWork_.nP); break; } @@ -411,16 +419,16 @@ void MahiFit::nnls() const { // no realizable optimization here (because it autovectorizes!) double minratio = std::numeric_limits::max(); for (unsigned int ipulse=0; ipulse( new ROOT::Math::Functor(psfPtr_.get(),&FitterFuncs::PulseShapeFunctor::singlePulseShapeFunc, 3) ); - } @@ -487,6 +493,7 @@ void MahiFit::nnlsUnconstrainParameter(Index idxp) const { nnlsWork_.aTaMat.col(nnlsWork_.nP).swap(nnlsWork_.aTaMat.col(idxp)); nnlsWork_.aTaMat.row(nnlsWork_.nP).swap(nnlsWork_.aTaMat.row(idxp)); nnlsWork_.pulseMat.col(nnlsWork_.nP).swap(nnlsWork_.pulseMat.col(idxp)); + nnlsWork_.pulseDerivMat.col(nnlsWork_.nP).swap(nnlsWork_.pulseDerivMat.col(idxp)); std::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP),nnlsWork_.aTbVec.coeffRef(idxp)); std::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP),nnlsWork_.ampVec.coeffRef(idxp)); std::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP),nnlsWork_.bxs.coeffRef(idxp)); @@ -497,6 +504,7 @@ void MahiFit::nnlsConstrainParameter(Index minratioidx) const { nnlsWork_.aTaMat.col(nnlsWork_.nP-1).swap(nnlsWork_.aTaMat.col(minratioidx)); nnlsWork_.aTaMat.row(nnlsWork_.nP-1).swap(nnlsWork_.aTaMat.row(minratioidx)); nnlsWork_.pulseMat.col(nnlsWork_.nP-1).swap(nnlsWork_.pulseMat.col(minratioidx)); + nnlsWork_.pulseDerivMat.col(nnlsWork_.nP-1).swap(nnlsWork_.pulseDerivMat.col(minratioidx)); std::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP-1),nnlsWork_.aTbVec.coeffRef(minratioidx)); std::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP-1),nnlsWork_.ampVec.coeffRef(minratioidx)); std::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP-1),nnlsWork_.bxs.coeffRef(minratioidx)); @@ -656,16 +664,8 @@ void MahiFit::resetWorkspace() const { nnlsWork_.maxoffset=0; nnlsWork_.dt=0; - std::fill(std::begin(nnlsWork_.pulseN), std::end(nnlsWork_.pulseN), 0); - std::fill(std::begin(nnlsWork_.pulseM), std::end(nnlsWork_.pulseM), 0); - std::fill(std::begin(nnlsWork_.pulseP), std::end(nnlsWork_.pulseP), 0); - nnlsWork_.amplitudes.setZero(); - nnlsWork_.invCovMat.setZero(); nnlsWork_.noiseTerms.setZero(); nnlsWork_.pedConstraint.setZero(); - nnlsWork_.residuals.setZero(); - - } diff --git a/RecoLocalCalo/HcalRecAlgos/src/PulseShapeFunctor.cc b/RecoLocalCalo/HcalRecAlgos/src/PulseShapeFunctor.cc index f3c72d710a5dc..6b6c14e0e4d89 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/PulseShapeFunctor.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/PulseShapeFunctor.cc @@ -103,7 +103,17 @@ namespace FitterFuncs{ PulseShapeFunctor::~PulseShapeFunctor() { } - double PulseShapeFunctor::EvalPulse(const double *pars, const unsigned nPars) { + void PulseShapeFunctor::EvalPulse(const double *pars) { + + int time = (pars[0]+timeShift_-timeMean_)*HcalConst::invertnsPerBx; + funcShape(pulse_shape_, pars[0],pars[1],psFit_slew[time]); + + return; + + } + + + double PulseShapeFunctor::EvalPulseM2(const double *pars, const unsigned nPars) { unsigned i =0, j=0; @@ -176,16 +186,20 @@ namespace FitterFuncs{ return chisq; } + void PulseShapeFunctor::singlePulseShapeFuncMahi( const double *x ) { + return EvalPulse(x); + } + double PulseShapeFunctor::singlePulseShapeFunc( const double *x ) { - return EvalPulse(x,3); + return EvalPulseM2(x,3); } double PulseShapeFunctor::doublePulseShapeFunc( const double *x ) { - return EvalPulse(x,5); + return EvalPulseM2(x,5); } double PulseShapeFunctor::triplePulseShapeFunc( const double *x ) { - return EvalPulse(x,7); + return EvalPulseM2(x,7); } }