MahiFit Class Reference

#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	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_ =nullptr
const HcalTimeSlew *	hcalTimeSlewDelay_ =nullptr

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 constexpr int	pedestalBX_ = 100

Detailed Description

Definition at line 84 of file MahiFit.h.

Member Typedef Documentation

typedef BXVector::Index MahiFit::Index

Definition at line 107 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 88 of file MahiFit.h.

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

{ };

Member Function Documentation

double MahiFit::calculateArrivalTime ( ) const

private

Definition at line 327 of file MahiFit.cc.

References MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::bxs, MahiNnlsWorkspace::fullTSOffset, nnlsWork_, MahiNnlsWorkspace::nPulseTot, PFRecoTauDiscriminationByIsolation_cfi::offset, pedestalBX_, MahiNnlsWorkspace::pulseDerivArray, MahiNnlsWorkspace::pulseDerivMat, MahiNnlsWorkspace::residuals, and MahiNnlsWorkspace::tsSize.

Referenced by doFit().

                                           {

  nnlsWork_.pulseDerivMat.resize(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_.fullTSOffset-offset, nnlsWork_.tsSize);
    }
  }

  PulseVector solution = nnlsWork_.pulseDerivMat.colPivHouseholderQr().solve(nnlsWork_.residuals);
  return solution.coeff(itIndex)/nnlsWork_.ampVec.coeff(itIndex);

}

double MahiFit::calculateChiSq ( ) const

private

Definition at line 463 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.

References activeBXs_, MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, MahiNnlsWorkspace::bxOffset, bxOffsetConf_, MahiNnlsWorkspace::bxs, bxSizeConf_, calculateArrivalTime(), dynamicPed_, MahiNnlsWorkspace::errVec, MahiNnlsWorkspace::fullTSOffset, MaxFSVSize, minimize(), nnlsWork_, MahiNnlsWorkspace::nPulseTot, PFRecoTauDiscriminationByIsolation_cfi::offset, pedestalBX_, MahiNnlsWorkspace::pulseCovArray, MahiNnlsWorkspace::pulseDerivArray, MahiNnlsWorkspace::pulseMat, MahiNnlsWorkspace::pulseShapeArray, MahiNnlsWorkspace::residuals, BXVector< T >::resize(), MahiNnlsWorkspace::tsOffset, MahiNnlsWorkspace::tsSize, and updatePulseShape().

Referenced by phase1Apply().

                                                                     {

  unsigned int bxSize=1;

  if (nbx==1) {
    nnlsWork_.bxOffset = 0;
  }
  else {
    bxSize = bxSizeConf_;
    nnlsWork_.bxOffset = bxOffsetConf_;
  }

  nnlsWork_.bxs.resize(bxSize);

  if (nbx==1) {
    nnlsWork_.bxs.coeffRef(0) = 0;
  }
  else {
    for (unsigned int iBX=0; iBX<bxSize; ++iBX) {
      nnlsWork_.bxs.coeffRef(iBX) = activeBXs_[iBX];
    }
  }

  nnlsWork_.nPulseTot = bxSize;

  if (dynamicPed_) {
    nnlsWork_.nPulseTot++;
    nnlsWork_.bxs.resize(nnlsWork_.nPulseTot);
    nnlsWork_.bxs[nnlsWork_.nPulseTot-1] = pedestalBX_;
  }

  nnlsWork_.pulseMat.resize(nnlsWork_.tsSize,nnlsWork_.nPulseTot);
  nnlsWork_.ampVec = PulseVector::Zero(nnlsWork_.nPulseTot);
  nnlsWork_.errVec = PulseVector::Zero(nnlsWork_.nPulseTot);

  int offset=0;
  for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
    offset=nnlsWork_.bxs.coeff(iBX);

    nnlsWork_.pulseShapeArray[iBX] = FullSampleVector::Zero(MaxFSVSize);
    nnlsWork_.pulseDerivArray[iBX] = FullSampleVector::Zero(MaxFSVSize);
    nnlsWork_.pulseCovArray[iBX]   = FullSampleMatrix::Constant(0);

    if (offset==pedestalBX_) {
      nnlsWork_.ampVec.coeffRef(iBX) = 0;
    }
    else {

      updatePulseShape(nnlsWork_.amplitudes.coeff(nnlsWork_.tsOffset + offset), 
               nnlsWork_.pulseShapeArray[iBX], 
               nnlsWork_.pulseDerivArray[iBX],
               nnlsWork_.pulseCovArray[iBX]);
      
      nnlsWork_.ampVec.coeffRef(iBX)=0;

      nnlsWork_.pulseMat.col(iBX) = nnlsWork_.pulseShapeArray[iBX].segment(nnlsWork_.fullTSOffset - offset, nnlsWork_.tsSize);

    }
  }

  nnlsWork_.pulseMat.col(nnlsWork_.nPulseTot-1) = SampleVector::Ones(nnlsWork_.tsSize);

  nnlsWork_.aTaMat.resize(nnlsWork_.nPulseTot, nnlsWork_.nPulseTot);
  nnlsWork_.aTbVec.resize(nnlsWork_.nPulseTot);

  double chiSq = minimize(); 

  nnlsWork_.residuals = nnlsWork_.pulseMat*nnlsWork_.ampVec - nnlsWork_.amplitudes;

  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
    double arrivalTime = calculateArrivalTime();
    correctedOutput.at(1) = arrivalTime; //time
    correctedOutput.at(2) = chiSq; //chi2

  }
  
}

double MahiFit::minimize ( ) const

private

Definition at line 217 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 351 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::ampvecpermtest, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::bxs, MahiNnlsWorkspace::covDecomp, MahiNnlsWorkspace::fullTSOffset, mps_fire::i, MahiNnlsWorkspace::invcovp, SiStripPI::max, min(), nMaxItersNNLS_, nnlsConstrainParameter(), nnlsThresh_, nnlsUnconstrainParameter(), nnlsWork_, MahiNnlsWorkspace::nP, MahiNnlsWorkspace::nPulseTot, PFRecoTauDiscriminationByIsolation_cfi::offset, pedestalBX_, MahiNnlsWorkspace::pulseMat, MahiNnlsWorkspace::pulseShapeArray, particleFlowDisplacedVertex_cfi::ratio, solveSubmatrix(), electronIdCutBased_cfi::threshold, MahiNnlsWorkspace::tsSize, and MahiNnlsWorkspace::updateWork.

Referenced by minimize().

                         {
  const unsigned int npulse = nnlsWork_.nPulseTot;
  
  for (unsigned int iBX=0; iBX<npulse; ++iBX) {
    int offset=nnlsWork_.bxs.coeff(iBX);
    if (offset==pedestalBX_) {
      nnlsWork_.pulseMat.col(iBX) = SampleVector::Ones(nnlsWork_.tsSize);
    }
    else {
      nnlsWork_.pulseMat.col(iBX) = nnlsWork_.pulseShapeArray.at(offset+nnlsWork_.bxOffset).segment(nnlsWork_.fullTSOffset-offset, nnlsWork_.tsSize);
    }
  }

  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 503 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 493 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 451 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 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 = pedVal*SampleMatrix::Ones(nnlsWork_.tsSize, nnlsWork_.tsSize);
  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::resetPulseShapeTemplate

(

const HcalPulseShapes::Shape &

ps

)

Definition at line 481 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 583 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::ampvecpermtest, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, begin, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::bxs, MahiNnlsWorkspace::covDecompLinv, MahiNnlsWorkspace::dt, end, MahiNnlsWorkspace::errVec, lumiContext::fill, MahiNnlsWorkspace::fullTSOffset, MahiNnlsWorkspace::invCovMat, MahiNnlsWorkspace::invcovp, nnlsWork_, MahiNnlsWorkspace::noiseTerms, MahiNnlsWorkspace::nPulseTot, MahiNnlsWorkspace::pedConstraint, MahiNnlsWorkspace::pulseCovArray, MahiNnlsWorkspace::pulseDerivArray, MahiNnlsWorkspace::pulseDerivMat, MahiNnlsWorkspace::pulseM, MahiNnlsWorkspace::pulseMat, MahiNnlsWorkspace::pulseN, MahiNnlsWorkspace::pulseP, MahiNnlsWorkspace::pulseShapeArray, MahiNnlsWorkspace::residuals, MahiNnlsWorkspace::topleft_work, MahiNnlsWorkspace::tsOffset, MahiNnlsWorkspace::tsSize, and MahiNnlsWorkspace::updateWork.

Referenced by phase1Apply().

                                   {

  nnlsWork_.nPulseTot=0;
  nnlsWork_.tsSize=0;
  nnlsWork_.tsOffset=0;
  nnlsWork_.fullTSOffset=0;
  nnlsWork_.bxOffset=0;
  nnlsWork_.dt=0;

  std::fill(std::begin(nnlsWork_.pulseCovArray), std::end(nnlsWork_.pulseCovArray), FullSampleMatrix::Zero());
  std::fill(std::begin(nnlsWork_.pulseShapeArray), std::end(nnlsWork_.pulseShapeArray), FullSampleVector::Zero());
  std::fill(std::begin(nnlsWork_.pulseDerivArray), std::end(nnlsWork_.pulseDerivArray), FullSampleVector::Zero());

  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_.bxs.setZero();
  nnlsWork_.invCovMat.setZero();
  nnlsWork_.noiseTerms.setZero();
  nnlsWork_.pedConstraint.setZero();
  nnlsWork_.pulseMat.setZero();
  nnlsWork_.pulseDerivMat.setZero();
  nnlsWork_.residuals.setZero();
  nnlsWork_.ampVec.setZero();
  nnlsWork_.errVec.setZero();
  nnlsWork_.ampvecpermtest.setZero();
  nnlsWork_.invcovp.setZero();
  nnlsWork_.aTaMat.setZero();
  nnlsWork_.aTbVec.setZero();
  nnlsWork_.updateWork.setZero();
  nnlsWork_.covDecompLinv.setZero();
  nnlsWork_.topleft_work.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 468 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 514 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 305 of file MahiFit.cc.

References MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::bxs, MahiNnlsWorkspace::covDecomp, MahiNnlsWorkspace::fullTSOffset, MahiNnlsWorkspace::invCovMat, 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_.fullTSOffset-offset, nnlsWork_.fullTSOffset-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 250 of file MahiFit.cc.

References applyTimeSlew_, HcalTimeSlew::delay(), delta, MahiNnlsWorkspace::dt, MahiNnlsWorkspace::fullTSOffset, hcalTimeSlewDelay_, 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 =nnlsWork_. tsOffset == 3 ? 1 : 0;

  for (unsigned int iTS=nnlsWork_.fullTSOffset; iTS<nnlsWork_.fullTSOffset + nnlsWork_.tsSize; ++iTS) {

    pulseShape.coeffRef(iTS) = nnlsWork_.pulseN[iTS-nnlsWork_.fullTSOffset+delta];
    pulseDeriv.coeffRef(iTS) = 0.5*(nnlsWork_.pulseM[iTS-nnlsWork_.fullTSOffset+delta]+nnlsWork_.pulseP[iTS-nnlsWork_.fullTSOffset+delta])/(2*nnlsWork_.dt);

    nnlsWork_.pulseM[iTS-nnlsWork_.fullTSOffset] -= nnlsWork_.pulseN[iTS-nnlsWork_.fullTSOffset];
    nnlsWork_.pulseP[iTS-nnlsWork_.fullTSOffset] -= nnlsWork_.pulseN[iTS-nnlsWork_.fullTSOffset];
  }

  for (unsigned int iTS=nnlsWork_.fullTSOffset; iTS<nnlsWork_.fullTSOffset+nnlsWork_.tsSize; ++iTS) {
    for (unsigned int jTS=nnlsWork_.fullTSOffset; jTS<iTS+1; ++jTS) {
      
      double tmp = 0.5*( nnlsWork_.pulseP[iTS-nnlsWork_.fullTSOffset+delta]*nnlsWork_.pulseP[jTS-nnlsWork_.fullTSOffset+delta] +
             nnlsWork_.pulseM[iTS-nnlsWork_.fullTSOffset+delta]*nnlsWork_.pulseM[jTS-nnlsWork_.fullTSOffset+delta] );
      
      pulseCov(iTS,jTS) += tmp;
      pulseCov(jTS,iTS) += tmp;
      
    }
  }
  
}

Member Data Documentation

std::vector<int> MahiFit::activeBXs_

private

Definition at line 151 of file MahiFit.h.

Referenced by doFit(), and setParameters().

bool MahiFit::applyTimeSlew_

private

Definition at line 143 of file MahiFit.h.

Referenced by setParameters(), and updatePulseShape().

int MahiFit::bxOffsetConf_

private

Definition at line 160 of file MahiFit.h.

Referenced by doFit(), and setParameters().

unsigned int MahiFit::bxSizeConf_

private

Definition at line 159 of file MahiFit.h.

Referenced by doFit(), and setParameters().

float MahiFit::chiSqSwitch_

private

Definition at line 141 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

int MahiFit::cntsetPulseShape_

private

Definition at line 163 of file MahiFit.h.

Referenced by resetPulseShapeTemplate().

const HcalPulseShapes::Shape* MahiFit::currentPulseShape_ =nullptr

Definition at line 108 of file MahiFit.h.

Referenced by setPulseShapeTemplate().

float MahiFit::deltaChiSqThresh_

private

Definition at line 156 of file MahiFit.h.

Referenced by minimize(), and setParameters().

bool MahiFit::dynamicPed_

private

Definition at line 139 of file MahiFit.h.

Referenced by doFit(), and setParameters().

const unsigned int MahiFit::fullTSofInterest_

private

Definition at line 134 of file MahiFit.h.

Referenced by phase1Apply().

const unsigned int MahiFit::fullTSSize_

private

Definition at line 133 of file MahiFit.h.

const HcalTimeSlew* MahiFit::hcalTimeSlewDelay_ =nullptr

Definition at line 109 of file MahiFit.h.

Referenced by setPulseShapeTemplate(), and updatePulseShape().

float MahiFit::meanTime_

private

Definition at line 147 of file MahiFit.h.

Referenced by setParameters(), and updatePulseShape().

int MahiFit::nMaxItersMin_

private

Definition at line 153 of file MahiFit.h.

Referenced by minimize(), and setParameters().

int MahiFit::nMaxItersNNLS_

private

Definition at line 154 of file MahiFit.h.

Referenced by nnls(), and setParameters().

float MahiFit::nnlsThresh_

private

Definition at line 157 of file MahiFit.h.

Referenced by nnls(), and setParameters().

MahiNnlsWorkspace MahiFit::nnlsWork_

mutableprivate

Definition at line 130 of file MahiFit.h.

Referenced by calculateArrivalTime(), calculateChiSq(), doFit(), minimize(), nnls(), nnlsConstrainParameter(), nnlsUnconstrainParameter(), onePulseMinimize(), phase1Apply(), resetWorkspace(), updateCov(), and updatePulseShape().

constexpr int MahiFit::pedestalBX_ = 100

staticprivate

Definition at line 136 of file MahiFit.h.

Referenced by calculateArrivalTime(), doFit(), nnls(), and updateCov().

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

private

Definition at line 165 of file MahiFit.h.

Referenced by resetPulseShapeTemplate().

std::unique_ptr<FitterFuncs::PulseShapeFunctor> MahiFit::psfPtr_

private

Definition at line 164 of file MahiFit.h.

Referenced by resetPulseShapeTemplate(), and updatePulseShape().

HcalTimeSlew::BiasSetting MahiFit::slewFlavor_

private

Definition at line 144 of file MahiFit.h.

Referenced by setParameters(), setPulseShapeTemplate(), and updatePulseShape().

float MahiFit::timeSigmaHPD_

private

Definition at line 148 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

float MahiFit::timeSigmaSiPM_

private

Definition at line 149 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

float MahiFit::ts4Thresh_

private

Definition at line 140 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

double MahiFit::tsDelay1GeV_ =0

private

Definition at line 145 of file MahiFit.h.

Referenced by setPulseShapeTemplate(), and updatePulseShape().