#include <MahiFit.h>

Public Types
typedef BXVector::Index	Index

Public Member Functions
void	doFit (std::array< float, 3 > &correctedOutput, const int nbx) const

	MahiFit ()

void	phase1Apply (const HBHEChannelInfo &channelData, float &reconstructedEnergy, float &reconstructedTime, bool &useTriple, float &chi2) const

void	phase1Debug (const HBHEChannelInfo &channelData, MahiDebugInfo &mdi) const

void	resetPulseShapeTemplate (const HcalPulseShapes::Shape &ps)

void	setParameters (bool iDynamicPed, double iTS4Thresh, double chiSqSwitch, bool iApplyTimeSlew, HcalTimeSlew::BiasSetting slewFlavor, double iMeanTime, double iTimeSigmaHPD, double iTimeSigmaSiPM, const std::vector< int > &iActiveBXs, int iNMaxItersMin, int iNMaxItersNNLS, double iDeltaChiSqThresh, double iNnlsThresh)

void	setPulseShapeTemplate (const HcalPulseShapes::Shape &ps, const HcalTimeSlew *hcalTimeSlewDelay)

	~MahiFit ()

Public Attributes
const HcalPulseShapes::Shape *	currentPulseShape_ = 0

const HcalTimeSlew *	hcalTimeSlewDelay_ = 0

Private Member Functions
double	calculateArrivalTime () const

double	calculateChiSq () const

double	minimize () const

void	nnls () const

void	nnlsConstrainParameter (Index minratioidx) const

void	nnlsUnconstrainParameter (Index idxp) const

void	onePulseMinimize () const

void	resetWorkspace () const

void	solveSubmatrix (PulseMatrix &mat, PulseVector &invec, PulseVector &outvec, unsigned nP) const

void	updateCov () const

void	updatePulseShape (double itQ, FullSampleVector &pulseShape, FullSampleVector &pulseDeriv, FullSampleMatrix &pulseCov) const

Private Attributes
std::vector< int >	activeBXs_

bool	applyTimeSlew_

int	bxOffsetConf_

unsigned int	bxSizeConf_

float	chiSqSwitch_

int	cntsetPulseShape_

float	deltaChiSqThresh_

bool	dynamicPed_

const unsigned int	fullTSofInterest_

const unsigned int	fullTSSize_

float	meanTime_

int	nMaxItersMin_

int	nMaxItersNNLS_

float	nnlsThresh_

MahiNnlsWorkspace	nnlsWork_

std::unique_ptr< ROOT::Math::Functor >	pfunctor_

std::unique_ptr< FitterFuncs::PulseShapeFunctor >	psfPtr_

HcalTimeSlew::BiasSetting	slewFlavor_

float	timeSigmaHPD_

float	timeSigmaSiPM_

float	ts4Thresh_

double	tsDelay1GeV_ =0

Static Private Attributes
static int	pedestalBX_ = 100

static float	timeLimit_ = 12.5

Detailed Description

Definition at line 118 of file MahiFit.h.

Member Typedef Documentation

typedef BXVector::Index MahiFit::Index

Definition at line 144 of file MahiFit.h.

Constructor & Destructor Documentation

MahiFit::MahiFit ( )

Definition at line 4 of file MahiFit.cc.

                  :
   fullTSSize_(19), 
   fullTSofInterest_(8)
 {}

MahiFit::~MahiFit ( )

inline

Definition at line 122 of file MahiFit.h.

References vertices_cff::chi2, and HBHEMahiParameters_cfi::chiSqSwitch.

122 { };

Member Function Documentation

double MahiFit::calculateArrivalTime ( ) const

private

Definition at line 323 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::bxs, MahiNnlsWorkspace::maxoffset, nnlsWork_, MahiNnlsWorkspace::nPulseTot, PFRecoTauDiscriminationByIsolation_cfi::offset, pedestalBX_, MahiNnlsWorkspace::pulseDerivArray, MahiNnlsWorkspace::pulseDerivMat, MahiNnlsWorkspace::pulseMat, MahiNnlsWorkspace::residuals, lumiQTWidget::t, timeLimit_, and MahiNnlsWorkspace::tsSize.

Referenced by doFit().

                                            {
 
   nnlsWork_.residuals = nnlsWork_.pulseMat*nnlsWork_.ampVec - nnlsWork_.amplitudes;
 
   nnlsWork_.pulseDerivMat.setZero(nnlsWork_.tsSize,nnlsWork_.nPulseTot);
 
   int itIndex=0;
 
   for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
     int offset=nnlsWork_.bxs.coeff(iBX);
     if (offset==0) itIndex=iBX;
 
     if (offset==pedestalBX_) {
       nnlsWork_.pulseDerivMat.col(iBX) = SampleVector::Zero(nnlsWork_.tsSize);
     }
     else {
       nnlsWork_.pulseDerivMat.col(iBX) = nnlsWork_.pulseDerivArray.at(offset+nnlsWork_.bxOffset).segment(nnlsWork_.maxoffset-offset, nnlsWork_.tsSize);
     }
   }
 
   PulseVector solution = nnlsWork_.pulseDerivMat.colPivHouseholderQr().solve(nnlsWork_.residuals);
   float t = solution.coeff(itIndex)/nnlsWork_.ampVec.coeff(itIndex);
   t = (t>timeLimit_) ?  timeLimit_ : 
     ((t<-timeLimit_) ? -timeLimit_ : t);
 
   return t;
 
 }

double MahiFit::calculateChiSq ( ) const

private

Definition at line 456 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::covDecomp, nnlsWork_, and MahiNnlsWorkspace::pulseMat.

Referenced by minimize().

                                      {
   
   return (nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat*nnlsWork_.ampVec - nnlsWork_.amplitudes)).squaredNorm();
 }

void MahiFit::doFit	(	std::array< float, 3 > &	correctedOutput,
		const int	nbx
	)		const

Definition at line 128 of file MahiFit.cc.

Referenced by phase1Apply().

                                                                      {
 
   unsigned int bxSize=1;
 
   if (nbx==1) {
     nnlsWork_.bxOffset = 0;
   }
   else {
     bxSize = bxSizeConf_;
     nnlsWork_.bxOffset = static_cast<int>(nnlsWork_.tsOffset) >= bxOffsetConf_ ? bxOffsetConf_ : nnlsWork_.tsOffset;
   }
 
   nnlsWork_.nPulseTot = bxSize;
 
   if (dynamicPed_) nnlsWork_.nPulseTot++;
   nnlsWork_.bxs.setZero(nnlsWork_.nPulseTot);
 
   if (nbx==1) {
     nnlsWork_.bxs.coeffRef(0) = 0;
   }
   else {
     for (unsigned int iBX=0; iBX<bxSize; ++iBX) {
       nnlsWork_.bxs.coeffRef(iBX) = activeBXs_[iBX] - ((static_cast<int>(nnlsWork_.tsOffset) + activeBXs_[0]) >= 0 ? 0 : (nnlsWork_.tsOffset + activeBXs_[0]));
     }
   }
 
   nnlsWork_.maxoffset = nnlsWork_.bxs.coeffRef(bxSize-1);
   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);
 
   int offset=0;
   for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
     offset=nnlsWork_.bxs.coeff(iBX);
 
     nnlsWork_.pulseShapeArray[iBX].setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
     nnlsWork_.pulseDerivArray[iBX].setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
     nnlsWork_.pulseCovArray[iBX].setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset, nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
 
 
     if (offset==pedestalBX_) {
       nnlsWork_.pulseMat.col(iBX) = SampleVector::Ones(nnlsWork_.tsSize);
     }
     else {
       updatePulseShape(nnlsWork_.amplitudes.coeff(nnlsWork_.tsOffset + offset), 
                nnlsWork_.pulseShapeArray[iBX], 
                nnlsWork_.pulseDerivArray[iBX],
                nnlsWork_.pulseCovArray[iBX]);
       
       nnlsWork_.pulseMat.col(iBX) = nnlsWork_.pulseShapeArray[iBX].segment(nnlsWork_.maxoffset - offset, nnlsWork_.tsSize);
     }
   }
 
   nnlsWork_.aTaMat.setZero(nnlsWork_.nPulseTot, nnlsWork_.nPulseTot);
   nnlsWork_.aTbVec.setZero(nnlsWork_.nPulseTot);
   nnlsWork_.updateWork.setZero(nnlsWork_.nPulseTot);
 
   double chiSq = minimize(); 
 
   bool foundintime = false;
   unsigned int ipulseintime = 0;
 
   for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
     if (nnlsWork_.bxs.coeff(iBX)==0) {
       ipulseintime = iBX;
       foundintime = true;
     }
   }
 
   if (foundintime) {
     correctedOutput.at(0) = nnlsWork_.ampVec.coeff(ipulseintime); //charge
     if (correctedOutput.at(0)!=0) {
     double arrivalTime = calculateArrivalTime();
     correctedOutput.at(1) = arrivalTime; //time
     }
     else correctedOutput.at(1) = -9999;//time
 
     correctedOutput.at(2) = chiSq; //chi2
 
   }
   
 }

double MahiFit::minimize ( ) const

private

Definition at line 214 of file MahiFit.cc.

References funct::abs(), calculateChiSq(), deltaChiSqThresh_, nMaxItersMin_, nnls(), nnlsWork_, MahiNnlsWorkspace::nPulseTot, onePulseMinimize(), and updateCov().

Referenced by doFit().

                                {
 
   double oldChiSq=9999;
   double chiSq=oldChiSq;
 
   for( int iter=1; iter<nMaxItersMin_ ; ++iter) {
 
     updateCov();
 
     if (nnlsWork_.nPulseTot>1) {
       nnls();
     }
     else {
       onePulseMinimize();
     }
 
     double newChiSq=calculateChiSq();
     double deltaChiSq = newChiSq - chiSq;
 
     if (newChiSq==oldChiSq && newChiSq<chiSq) {
       break;
     }
     oldChiSq=chiSq;
     chiSq = newChiSq;
 
     if (std::abs(deltaChiSq)<deltaChiSqThresh_) break;
 
   }
 
   return chiSq;
 
 }

void MahiFit::nnls ( ) const

private

Definition at line 353 of file MahiFit.cc.

Referenced by minimize().

                          {
   const unsigned int npulse = nnlsWork_.nPulseTot;
 
   nnlsWork_.invcovp = nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat);
   nnlsWork_.aTaMat = nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.invcovp);
   nnlsWork_.aTbVec = nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes));
   
   int iter = 0;
   Index idxwmax = 0;
   double wmax = 0.0;
   double threshold = nnlsThresh_;
 
   nnlsWork_.nP=0;
   
   while (true) {    
     if (iter>0 || nnlsWork_.nP==0) {
       if ( nnlsWork_.nP==std::min(npulse, nnlsWork_.tsSize)) break;
       
       const unsigned int nActive = npulse - nnlsWork_.nP;
       nnlsWork_.updateWork = nnlsWork_.aTbVec - nnlsWork_.aTaMat*nnlsWork_.ampVec;
       
       Index idxwmaxprev = idxwmax;
       double wmaxprev = wmax;
       wmax = nnlsWork_.updateWork.tail(nActive).maxCoeff(&idxwmax);
       
       if (wmax<threshold || (idxwmax==idxwmaxprev && wmax==wmaxprev)) {
     break;
       }
       
       if (iter>=nMaxItersNNLS_) {
     break;
       }
 
       //unconstrain parameter
       Index idxp = nnlsWork_.nP + idxwmax;
       nnlsUnconstrainParameter(idxp);
 
     }
 
     while (true) {
       if (nnlsWork_.nP==0) break;     
 
       nnlsWork_.ampvecpermtest = nnlsWork_.ampVec;
       
       solveSubmatrix(nnlsWork_.aTaMat,nnlsWork_.aTbVec,nnlsWork_.ampvecpermtest,nnlsWork_.nP);
 
       //check solution
       bool positive = true;
       for (unsigned int i = 0; i < nnlsWork_.nP; ++i)
         positive &= (nnlsWork_.ampvecpermtest(i) > 0);
       if (positive) {
         nnlsWork_.ampVec.head(nnlsWork_.nP) = nnlsWork_.ampvecpermtest.head(nnlsWork_.nP);
         break;
       } 
 
       //update parameter vector
       Index minratioidx=0;
       
       // no realizable optimization here (because it autovectorizes!)
       double minratio = std::numeric_limits<double>::max();
       for (unsigned int ipulse=0; ipulse<nnlsWork_.nP; ++ipulse) {
     if (nnlsWork_.ampvecpermtest.coeff(ipulse)<=0.) {
       const double c_ampvec = nnlsWork_.ampVec.coeff(ipulse);
       const double ratio = c_ampvec/(c_ampvec-nnlsWork_.ampvecpermtest.coeff(ipulse));
       if (ratio<minratio) {
         minratio = ratio;
         minratioidx = ipulse;
       }
     }
       }
       nnlsWork_.ampVec.head(nnlsWork_.nP) += minratio*(nnlsWork_.ampvecpermtest.head(nnlsWork_.nP) - nnlsWork_.ampVec.head(nnlsWork_.nP));
       
       //avoid numerical problems with later ==0. check
       nnlsWork_.ampVec.coeffRef(minratioidx) = 0.;
       
       nnlsConstrainParameter(minratioidx);
     }
    
     ++iter;
 
     //adaptive convergence threshold to avoid infinite loops but still
     //ensure best value is used
     if (iter%10==0) {
       threshold *= 10.;
     }
     
   }
 
   
 }

void MahiFit::nnlsConstrainParameter ( Index minratioidx ) const

private

Definition at line 496 of file MahiFit.cc.

References MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, MahiNnlsWorkspace::bxs, nnlsWork_, MahiNnlsWorkspace::nP, MahiNnlsWorkspace::pulseMat, and std::swap().

Referenced by nnls().

                                                             {
   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));
   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));
   --nnlsWork_.nP;
 
 }

void MahiFit::nnlsUnconstrainParameter ( Index idxp ) const

private

Definition at line 486 of file MahiFit.cc.

References MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, MahiNnlsWorkspace::bxs, nnlsWork_, MahiNnlsWorkspace::nP, MahiNnlsWorkspace::pulseMat, and std::swap().

Referenced by nnls().

                                                        {
   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));
   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));
   ++nnlsWork_.nP;
 }

void MahiFit::onePulseMinimize ( ) const

private

Definition at line 444 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::covDecomp, MahiNnlsWorkspace::invcovp, SiStripPI::max, nnlsWork_, and MahiNnlsWorkspace::pulseMat.

Referenced by minimize().

                                      {
 
   nnlsWork_.invcovp = nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat);
 
   SingleMatrix aTamatval = nnlsWork_.invcovp.transpose()*nnlsWork_.invcovp;
   SingleVector aTbvecval = nnlsWork_.invcovp.transpose()*nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes);
 
   nnlsWork_.ampVec.coeffRef(0) = std::max(0., aTbvecval.coeff(0)/aTamatval.coeff(0));
 
 }

void MahiFit::phase1Apply	(	const HBHEChannelInfo &	channelData,
		float &	reconstructedEnergy,
		float &	reconstructedTime,
		bool &	useTriple,
		float &	chi2
	)		const

Definition at line 40 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, ALCARECOTkAlJpsiMuMu_cff::charge, chiSqSwitch_, doFit(), MahiNnlsWorkspace::dt, HBHEChannelInfo::fcByPE(), MahiNnlsWorkspace::fullTSOffset, fullTSofInterest_, HBHEChannelInfo::hasTimeInfo(), nnlsWork_, MahiNnlsWorkspace::noiseTerms, HBHEChannelInfo::nSamples(), MahiNnlsWorkspace::pedConstraint, resetWorkspace(), HBHEChannelInfo::soi(), mathSSE::sqrt(), timeSigmaHPD_, timeSigmaSiPM_, ts4Thresh_, HBHEChannelInfo::tsDFcPerADC(), HBHEChannelInfo::tsGain(), MahiNnlsWorkspace::tsOffset, HBHEChannelInfo::tsPedestal(), HBHEChannelInfo::tsPedestalWidth(), HBHEChannelInfo::tsRawCharge(), and MahiNnlsWorkspace::tsSize.

Referenced by phase1Debug(), and SimpleHBHEPhase1Algo::reconstruct().

                                  {
 
   assert(channelData.nSamples()==8||channelData.nSamples()==10);
 
   resetWorkspace();
 
   nnlsWork_.tsSize = channelData.nSamples();
   nnlsWork_.tsOffset = channelData.soi();
   nnlsWork_.fullTSOffset = fullTSofInterest_ - nnlsWork_.tsOffset;
 
   // 1 sigma time constraint
   if (channelData.hasTimeInfo()) nnlsWork_.dt=timeSigmaSiPM_;
   else nnlsWork_.dt=timeSigmaHPD_;
 
 
   //Average pedestal width (for covariance matrix constraint)
   float pedVal = 0.25*( channelData.tsPedestalWidth(0)*channelData.tsPedestalWidth(0)+
             channelData.tsPedestalWidth(1)*channelData.tsPedestalWidth(1)+
             channelData.tsPedestalWidth(2)*channelData.tsPedestalWidth(2)+
             channelData.tsPedestalWidth(3)*channelData.tsPedestalWidth(3) );
 
   nnlsWork_.pedConstraint.setConstant(nnlsWork_.tsSize, nnlsWork_.tsSize, pedVal);
   nnlsWork_.amplitudes.resize(nnlsWork_.tsSize);
   nnlsWork_.noiseTerms.resize(nnlsWork_.tsSize);
 
   std::array<float,3> reconstructedVals {{ 0.0, -9999, -9999 }};
   
   double tsTOT = 0, tstrig = 0; // in GeV
   for(unsigned int iTS=0; iTS<nnlsWork_.tsSize; ++iTS){
     double charge = channelData.tsRawCharge(iTS);
     double ped = channelData.tsPedestal(iTS);
 
     nnlsWork_.amplitudes.coeffRef(iTS) = charge - ped;
    
     //ADC granularity
     double noiseADC = (1./sqrt(12))*channelData.tsDFcPerADC(iTS);
 
     //Photostatistics
     double noisePhoto = 0;
     if ( (charge-ped)>channelData.tsPedestalWidth(iTS)) {
       noisePhoto = sqrt((charge-ped)*channelData.fcByPE());
     }
 
     //Electronic pedestal
     double pedWidth = channelData.tsPedestalWidth(iTS);
 
     //Total uncertainty from all sources
     nnlsWork_.noiseTerms.coeffRef(iTS) = noiseADC*noiseADC + noisePhoto*noisePhoto + pedWidth*pedWidth;
 
     tsTOT += (charge - ped)*channelData.tsGain(0);
     if( iTS==nnlsWork_.tsOffset ){
       tstrig += (charge - ped)*channelData.tsGain(0);
     }
   }
 
   if(tstrig >= ts4Thresh_ && tsTOT > 0) {
 
     useTriple=false;
 
     // only do pre-fit with 1 pulse if chiSq threshold is positive
     if (chiSqSwitch_>0) {
       doFit(reconstructedVals,1);
       if (reconstructedVals[2]>chiSqSwitch_) {
     doFit(reconstructedVals,0); //nbx=0 means use configured BXs
     useTriple=true;
       }
     }
     else {
       doFit(reconstructedVals,0);
       useTriple=true;
     }
   }
   else{
     reconstructedVals.at(0) = 0.; //energy
     reconstructedVals.at(1) = -9999.; //time
     reconstructedVals.at(2) = -9999.; //chi2
   }
   
   reconstructedEnergy = reconstructedVals[0]*channelData.tsGain(0);
   reconstructedTime = reconstructedVals[1];
   chi2 = reconstructedVals[2];
 
 }

void MahiFit::phase1Debug	(	const HBHEChannelInfo &	channelData,
		MahiDebugInfo &	mdi
	)		const

Definition at line 507 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiDebugInfo::arrivalTime, MahiNnlsWorkspace::bxs, ALCARECOTkAlJpsiMuMu_cff::charge, vertices_cff::chi2, MahiDebugInfo::chiSq, MahiDebugInfo::count, MahiNnlsWorkspace::dt, HBHEChannelInfo::fcByPE(), HBHEChannelInfo::hasTimeInfo(), MahiDebugInfo::inGain, MahiDebugInfo::inNoiseADC, MahiDebugInfo::inNoisePhoto, MahiDebugInfo::inPedAvg, MahiDebugInfo::inPedestal, MahiDebugInfo::inputTDC, MahiDebugInfo::inputTS, MahiDebugInfo::inTimeConst, MahiDebugInfo::itPulse, MahiDebugInfo::mahiEnergy, MahiDebugInfo::nEnergy, nnlsWork_, MahiNnlsWorkspace::noiseTerms, MahiDebugInfo::nPulse, MahiNnlsWorkspace::nPulseTot, MahiDebugInfo::nSamples, HBHEChannelInfo::nSamples(), MahiDebugInfo::pedEnergy, pedestalBX_, MahiDebugInfo::pEnergy, phase1Apply(), MahiDebugInfo::pPulse, MahiNnlsWorkspace::pulseMat, MahiDebugInfo::soi, HBHEChannelInfo::soi(), mathSSE::sqrt(), MahiDebugInfo::totalUCNoise, HBHEChannelInfo::tsDFcPerADC(), HBHEChannelInfo::tsGain(), HBHEChannelInfo::tsPedestal(), HBHEChannelInfo::tsPedestalWidth(), HBHEChannelInfo::tsRawCharge(), HBHEChannelInfo::tsRiseTime(), MahiNnlsWorkspace::tsSize, and MahiDebugInfo::use3.

                                         {
 
   float recoEnergy, recoTime, chi2;
   bool use3;
   phase1Apply(channelData, recoEnergy, recoTime, use3, chi2);
 
 
   mdi.nSamples    = channelData.nSamples();
   mdi.soi         = channelData.soi();
 
   mdi.use3        = use3;
 
   mdi.inTimeConst = nnlsWork_.dt;
   mdi.inPedAvg    = 0.25*( channelData.tsPedestalWidth(0)*channelData.tsPedestalWidth(0)+
                channelData.tsPedestalWidth(1)*channelData.tsPedestalWidth(1)+
                channelData.tsPedestalWidth(2)*channelData.tsPedestalWidth(2)+
                channelData.tsPedestalWidth(3)*channelData.tsPedestalWidth(3) );
   mdi.inGain      = channelData.tsGain(0);
 
   for (unsigned int iTS=0; iTS<channelData.nSamples(); ++iTS) {
 
     double charge = channelData.tsRawCharge(iTS);
     double ped = channelData.tsPedestal(iTS);
 
     mdi.inNoiseADC[iTS]  = (1./sqrt(12))*channelData.tsDFcPerADC(iTS);
 
     if ( (charge-ped)>channelData.tsPedestalWidth(iTS)) {
       mdi.inNoisePhoto[iTS] = sqrt((charge-ped)*channelData.fcByPE());
     }
     else { mdi.inNoisePhoto[iTS] = 0; }
 
     mdi.inPedestal[iTS]  = channelData.tsPedestalWidth(iTS);    
     mdi.totalUCNoise[iTS] = nnlsWork_.noiseTerms.coeffRef(iTS);
 
     if (channelData.hasTimeInfo()) {
       mdi.inputTDC[iTS] = channelData.tsRiseTime(iTS);
     }
     else { mdi.inputTDC[iTS]=-1; }
 
   }
 
   mdi.arrivalTime = recoTime;
   mdi.chiSq       = chi2;
 
   for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
     if (nnlsWork_.bxs.coeff(iBX)==0) {
       mdi.mahiEnergy=nnlsWork_.ampVec.coeff(iBX);
       for(unsigned int iTS=0; iTS<nnlsWork_.tsSize; ++iTS){
     mdi.count[iTS] = iTS;
     mdi.inputTS[iTS] = nnlsWork_.amplitudes.coeff(iTS);
     mdi.itPulse[iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
       }
     }
     else if (nnlsWork_.bxs.coeff(iBX)==pedestalBX_) {
       mdi.pedEnergy=nnlsWork_.ampVec.coeff(iBX);
     }
     else if (nnlsWork_.bxs.coeff(iBX)==-1) {
       mdi.pEnergy=nnlsWork_.ampVec.coeff(iBX);
       for(unsigned int iTS=0; iTS<nnlsWork_.tsSize; ++iTS){
         mdi.pPulse[iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
       }
     }
     else if (nnlsWork_.bxs.coeff(iBX)==1) {
       mdi.nEnergy=nnlsWork_.ampVec.coeff(iBX);
       for(unsigned int iTS=0; iTS<nnlsWork_.tsSize; ++iTS){
         mdi.nPulse[iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
       }
     }
   }  
 }

void MahiFit::resetPulseShapeTemplate ( const HcalPulseShapes::Shape & ps )

Definition at line 474 of file MahiFit.cc.

References cntsetPulseShape_, pfunctor_, psfPtr_, and FitterFuncs::PulseShapeFunctor::singlePulseShapeFunc().

Referenced by setPulseShapeTemplate().

                                                                     { 
   ++ cntsetPulseShape_;
 
   // only the pulse shape itself from PulseShapeFunctor is used for Mahi
   // the uncertainty terms calculated inside PulseShapeFunctor are used for Method 2 only
   psfPtr_.reset(new FitterFuncs::PulseShapeFunctor(ps,false,false,false,
                            1,0,0,10));
   pfunctor_ = std::unique_ptr<ROOT::Math::Functor>( new ROOT::Math::Functor(psfPtr_.get(),&FitterFuncs::PulseShapeFunctor::singlePulseShapeFunc, 3) );
 
 
 }

void MahiFit::resetWorkspace ( ) const

private

Definition at line 649 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, begin, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::dt, end, lumiContext::fill, MahiNnlsWorkspace::fullTSOffset, MahiNnlsWorkspace::invCovMat, MahiNnlsWorkspace::maxoffset, nnlsWork_, MahiNnlsWorkspace::noiseTerms, MahiNnlsWorkspace::nPulseTot, MahiNnlsWorkspace::pedConstraint, MahiNnlsWorkspace::pulseM, MahiNnlsWorkspace::pulseN, MahiNnlsWorkspace::pulseP, MahiNnlsWorkspace::residuals, MahiNnlsWorkspace::tsOffset, and MahiNnlsWorkspace::tsSize.

Referenced by phase1Apply().

                                    {
 
   nnlsWork_.nPulseTot=0;
   nnlsWork_.tsSize=0;
   nnlsWork_.tsOffset=0;
   nnlsWork_.fullTSOffset=0;
   nnlsWork_.bxOffset=0;
   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();
 
 
 
 }

void MahiFit::setParameters	(	bool	iDynamicPed,
		double	iTS4Thresh,
		double	chiSqSwitch,
		bool	iApplyTimeSlew,
		HcalTimeSlew::BiasSetting	slewFlavor,
		double	iMeanTime,
		double	iTimeSigmaHPD,
		double	iTimeSigmaSiPM,
		const std::vector< int > &	iActiveBXs,
		int	iNMaxItersMin,
		int	iNMaxItersNNLS,
		double	iDeltaChiSqThresh,
		double	iNnlsThresh
	)

Definition at line 9 of file MahiFit.cc.

References activeBXs_, applyTimeSlew_, bxOffsetConf_, bxSizeConf_, HBHEMahiParameters_cfi::chiSqSwitch, chiSqSwitch_, deltaChiSqThresh_, dynamicPed_, meanTime_, nMaxItersMin_, nMaxItersNNLS_, nnlsThresh_, slewFlavor_, timeSigmaHPD_, timeSigmaSiPM_, and ts4Thresh_.

                                                               {
 
   dynamicPed_    = iDynamicPed;
 
   ts4Thresh_     = iTS4Thresh;
   chiSqSwitch_   = chiSqSwitch;
 
   applyTimeSlew_ = iApplyTimeSlew;
   slewFlavor_    = slewFlavor;
 
   meanTime_      = iMeanTime;
   timeSigmaHPD_  = iTimeSigmaHPD;
   timeSigmaSiPM_ = iTimeSigmaSiPM;
 
   activeBXs_ = iActiveBXs;
 
   nMaxItersMin_  = iNMaxItersMin;
   nMaxItersNNLS_ = iNMaxItersNNLS;
 
   deltaChiSqThresh_ = iDeltaChiSqThresh;
   nnlsThresh_    = iNnlsThresh;
 
   bxOffsetConf_ = -(*std::min_element(activeBXs_.begin(), activeBXs_.end()));
   bxSizeConf_   = activeBXs_.size();
 
 }

void MahiFit::setPulseShapeTemplate	(	const HcalPulseShapes::Shape &	ps,
		const HcalTimeSlew *	hcalTimeSlewDelay
	)

Definition at line 461 of file MahiFit.cc.

References currentPulseShape_, HcalTimeSlew::delay(), hcalTimeSlewDelay_, resetPulseShapeTemplate(), slewFlavor_, and tsDelay1GeV_.

                                                                                                         {
 
   if (!(&ps == currentPulseShape_ ))
     {
 
       hcalTimeSlewDelay_ = hcalTimeSlewDelay;
       tsDelay1GeV_= hcalTimeSlewDelay->delay(1.0, slewFlavor_);
 
       resetPulseShapeTemplate(ps);
       currentPulseShape_ = &ps;
     }
 }

void MahiFit::solveSubmatrix	(	PulseMatrix &	mat,
		PulseVector &	invec,
		PulseVector &	outvec,
		unsigned	nP
	)		const

private

Definition at line 580 of file MahiFit.cc.

References Exception, and groupFilesInBlocks::temp.

Referenced by nnls().

                                                                                                          {
   using namespace Eigen;
   switch( nP ) { // pulse matrix is always square.
   case 10:
     {
       Matrix<double,10,10> temp = mat;
       outvec.head<10>() = temp.ldlt().solve(invec.head<10>());
     }
     break;
   case 9:
     {
       Matrix<double,9,9> temp = mat.topLeftCorner<9,9>();
       outvec.head<9>() = temp.ldlt().solve(invec.head<9>());
     }
     break;
   case 8:
     {
       Matrix<double,8,8> temp = mat.topLeftCorner<8,8>();
       outvec.head<8>() = temp.ldlt().solve(invec.head<8>());
     }
     break;
   case 7:
     {
       Matrix<double,7,7> temp = mat.topLeftCorner<7,7>();
       outvec.head<7>() = temp.ldlt().solve(invec.head<7>());
     }
     break;
   case 6:
     {
       Matrix<double,6,6> temp = mat.topLeftCorner<6,6>();
       outvec.head<6>() = temp.ldlt().solve(invec.head<6>());
     }
     break;
   case 5:
     {
       Matrix<double,5,5> temp = mat.topLeftCorner<5,5>();
       outvec.head<5>() = temp.ldlt().solve(invec.head<5>());
     }
     break;
   case 4:
     {
       Matrix<double,4,4> temp = mat.topLeftCorner<4,4>();
       outvec.head<4>() = temp.ldlt().solve(invec.head<4>());
     }
     break;
   case 3: 
     {
       Matrix<double,3,3> temp = mat.topLeftCorner<3,3>();
       outvec.head<3>() = temp.ldlt().solve(invec.head<3>());
     }
     break;
   case 2:
     {
       Matrix<double,2,2> temp = mat.topLeftCorner<2,2>();
       outvec.head<2>() = temp.ldlt().solve(invec.head<2>());
     }
     break;
   case 1:
     {
       Matrix<double,1,1> temp = mat.topLeftCorner<1,1>();
       outvec.head<1>() = temp.ldlt().solve(invec.head<1>());
     }
     break;
   default:
     throw cms::Exception("HcalMahiWeirdState") 
       << "Weird number of pulses encountered in Mahi, module is configured incorrectly!";
   }
 }

void MahiFit::updateCov ( ) const

private

Definition at line 301 of file MahiFit.cc.

References MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::bxs, MahiNnlsWorkspace::covDecomp, MahiNnlsWorkspace::invCovMat, MahiNnlsWorkspace::maxoffset, nnlsWork_, MahiNnlsWorkspace::noiseTerms, MahiNnlsWorkspace::nPulseTot, PFRecoTauDiscriminationByIsolation_cfi::offset, MahiNnlsWorkspace::pedConstraint, pedestalBX_, MahiNnlsWorkspace::pulseCovArray, and MahiNnlsWorkspace::tsSize.

Referenced by minimize().

                               {
 
   nnlsWork_.invCovMat.resize(nnlsWork_.tsSize, nnlsWork_.tsSize);
 
   nnlsWork_.invCovMat = nnlsWork_.noiseTerms.asDiagonal();
   nnlsWork_.invCovMat +=nnlsWork_.pedConstraint;
 
   for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
     if (nnlsWork_.ampVec.coeff(iBX)==0) continue;
     
     int offset=nnlsWork_.bxs.coeff(iBX);
 
     if (offset==pedestalBX_) continue;             
     else { 
       nnlsWork_.invCovMat += nnlsWork_.ampVec.coeff(iBX)*nnlsWork_.ampVec.coeff(iBX)
     *nnlsWork_.pulseCovArray.at(offset+nnlsWork_.bxOffset).block(nnlsWork_.maxoffset-offset, nnlsWork_.maxoffset-offset, nnlsWork_.tsSize, nnlsWork_.tsSize);
     }
   }
   
   nnlsWork_.covDecomp.compute(nnlsWork_.invCovMat);
 }

void MahiFit::updatePulseShape	(	double	itQ,
		FullSampleVector &	pulseShape,
		FullSampleVector &	pulseDeriv,
		FullSampleMatrix &	pulseCov
	)		const

private

Definition at line 247 of file MahiFit.cc.

References applyTimeSlew_, HcalTimeSlew::delay(), delta, MahiNnlsWorkspace::dt, hcalTimeSlewDelay_, MahiNnlsWorkspace::maxoffset, meanTime_, nnlsWork_, psfPtr_, MahiNnlsWorkspace::pulseM, MahiNnlsWorkspace::pulseN, MahiNnlsWorkspace::pulseP, slewFlavor_, cscNeutronWriter_cfi::t0, tmp, tsDelay1GeV_, MahiNnlsWorkspace::tsSize, and geometryCSVtoXML::xx.

Referenced by doFit().

                                                      {
   
   float t0=meanTime_;
 
   if(applyTimeSlew_) {
     if(itQ<=1.0) t0+=tsDelay1GeV_;
     else t0+=hcalTimeSlewDelay_->delay(itQ,slewFlavor_);
   }
 
   nnlsWork_.pulseN.fill(0);
   nnlsWork_.pulseM.fill(0);
   nnlsWork_.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);
 
   (*pfunctor_)(&xxm[0]);
   psfPtr_->getPulseShape(nnlsWork_.pulseM);
   
   (*pfunctor_)(&xxp[0]);
   psfPtr_->getPulseShape(nnlsWork_.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;
 
   for (unsigned int iTS=0; iTS<nnlsWork_.tsSize; ++iTS) {
 
     pulseShape.coeffRef(iTS+nnlsWork_.maxoffset) = nnlsWork_.pulseN[iTS+delta];
     pulseDeriv.coeffRef(iTS+nnlsWork_.maxoffset) = 0.5*(nnlsWork_.pulseM[iTS+delta]-nnlsWork_.pulseP[iTS+delta])/(2*nnlsWork_.dt);
 
     nnlsWork_.pulseM[iTS] -= nnlsWork_.pulseN[iTS];
     nnlsWork_.pulseP[iTS] -= nnlsWork_.pulseN[iTS];
   }
 
   for (unsigned int iTS=0; iTS<nnlsWork_.tsSize; ++iTS) {
     for (unsigned int jTS=0; jTS<iTS+1; ++jTS) {
       
       double tmp = 0.5*( nnlsWork_.pulseP[iTS+delta]*nnlsWork_.pulseP[jTS+delta] +
              nnlsWork_.pulseM[iTS+delta]*nnlsWork_.pulseM[jTS+delta] );
   
       pulseCov(jTS+nnlsWork_.maxoffset,iTS+nnlsWork_.maxoffset) += tmp;      
       if(iTS!=jTS) pulseCov(iTS+nnlsWork_.maxoffset,jTS+nnlsWork_.maxoffset) += tmp;
       
     }
   }
   
 }