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	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, bool iCalculateArrivalTime, 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
float	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_
bool	calculateArrivalTime_
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_
float	tsDelay1GeV_ =0.f

Static Private Attributes
static int	pedestalBX_ = 100
static float	timeLimit_ = 12.5

Detailed Description

Definition at line 98 of file MahiFit.h.

Member Typedef Documentation

typedef BXVector::Index MahiFit::Index

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

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

{ };

Member Function Documentation

float MahiFit::calculateArrivalTime ( ) const

private

Definition at line 335 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::bxs, nnlsWork_, MahiNnlsWorkspace::nPulseTot, PFRecoTauDiscriminationByIsolation_cfi::offset, MahiNnlsWorkspace::pulseDerivMat, MahiNnlsWorkspace::pulseMat, protons_cff::t, and timeLimit_.

Referenced by doFit().

                                          {

  int itIndex=0;

  for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
    int offset=nnlsWork_.bxs.coeff(iBX);
    if (offset==0) itIndex=iBX;
    nnlsWork_.pulseDerivMat.col(iBX) *= nnlsWork_.ampVec.coeff(iBX);

  }

  SampleVector residuals = nnlsWork_.pulseMat*nnlsWork_.ampVec - nnlsWork_.amplitudes;
  PulseVector solution = nnlsWork_.pulseDerivMat.colPivHouseholderQr().solve(residuals);
  float t = solution.coeff(itIndex);
  t = (t>timeLimit_) ?  timeLimit_ : 
    ((t<-timeLimit_) ? -timeLimit_ : t);

  return t;

}

double MahiFit::calculateChiSq ( ) const

private

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

References activeBXs_, MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::bxOffset, bxOffsetConf_, MahiNnlsWorkspace::bxs, bxSizeConf_, calculateArrivalTime(), calculateArrivalTime_, dynamicPed_, MahiNnlsWorkspace::maxoffset, minimize(), nnlsWork_, MahiNnlsWorkspace::nPulseTot, PFRecoTauDiscriminationByIsolation_cfi::offset, pedestalBX_, MahiNnlsWorkspace::pulseCovArray, MahiNnlsWorkspace::pulseDerivMat, MahiNnlsWorkspace::pulseMat, MahiNnlsWorkspace::tsOffset, MahiNnlsWorkspace::tsSize, and updatePulseShape().

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_.pulseDerivMat.setZero(nnlsWork_.tsSize,nnlsWork_.nPulseTot);

  FullSampleVector pulseShapeArray;
  FullSampleVector pulseDerivArray;

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

    pulseShapeArray.setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
    pulseDerivArray.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);
      nnlsWork_.pulseDerivMat.col(iBX) = SampleVector::Zero(nnlsWork_.tsSize);
    }
    else {
      updatePulseShape(nnlsWork_.amplitudes.coeff(nnlsWork_.tsOffset + offset), 
               pulseShapeArray,
               pulseDerivArray,
               nnlsWork_.pulseCovArray[iBX]);
      
      nnlsWork_.pulseMat.col(iBX) = pulseShapeArray.segment(nnlsWork_.maxoffset - offset, nnlsWork_.tsSize);
      nnlsWork_.pulseDerivMat.col(iBX) = pulseDerivArray.segment(nnlsWork_.maxoffset-offset, nnlsWork_.tsSize);
    }
  }

  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) {
      // fixME store the timeslew
    float arrivalTime = 0.;
    if(calculateArrivalTime_) 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(), MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, calculateChiSq(), deltaChiSqThresh_, MahiNnlsWorkspace::invcovp, nMaxItersMin_, nnls(), nnlsWork_, MahiNnlsWorkspace::nPulseTot, onePulseMinimize(), MahiNnlsWorkspace::tsSize, and updateCov().

Referenced by doFit().

                               {

  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;

  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 357 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, MahiNnlsWorkspace::covDecomp, mps_fire::i, MahiNnlsWorkspace::invcovp, SiStripPI::max, min(), nMaxItersNNLS_, nnlsConstrainParameter(), nnlsThresh_, nnlsUnconstrainParameter(), nnlsWork_, MahiNnlsWorkspace::nP, MahiNnlsWorkspace::nPulseTot, MahiNnlsWorkspace::pulseMat, particleFlowDisplacedVertex_cfi::ratio, solveSubmatrix(), electronIdCutBased_cfi::threshold, and MahiNnlsWorkspace::tsSize.

Referenced by minimize().

                         {
  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);
  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, nsamples)) break;
      
      const unsigned int nActive = npulse - nnlsWork_.nP;
      updateWork = nnlsWork_.aTbVec - nnlsWork_.aTaMat*nnlsWork_.ampVec;
      
      Index idxwmaxprev = idxwmax;
      double wmaxprev = wmax;
      wmax = 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;     

      ampvecpermtest = nnlsWork_.ampVec;
      
      solveSubmatrix(nnlsWork_.aTaMat,nnlsWork_.aTbVec,ampvecpermtest,nnlsWork_.nP);

      //check solution
      bool positive = true;
      for (unsigned int i = 0; i < nnlsWork_.nP; ++i)
        positive &= (ampvecpermtest(i) > 0);
      if (positive) {
        nnlsWork_.ampVec.head(nnlsWork_.nP) = 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 (ampvecpermtest.coeff(ipulse)<=0.) {
      const double c_ampvec = nnlsWork_.ampVec.coeff(ipulse);
      const double ratio = c_ampvec/(c_ampvec-ampvecpermtest.coeff(ipulse));
      if (ratio<minratio) {
        minratio = ratio;
        minratioidx = ipulse;
      }
    }
      }
      nnlsWork_.ampVec.head(nnlsWork_.nP) += minratio*(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 506 of file MahiFit.cc.

References MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, MahiNnlsWorkspace::bxs, nnlsWork_, MahiNnlsWorkspace::nP, MahiNnlsWorkspace::pulseDerivMat, 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));
  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));
  --nnlsWork_.nP;

}

void MahiFit::nnlsUnconstrainParameter ( Index  idxp ) const

private

Definition at line 495 of file MahiFit.cc.

References MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, MahiNnlsWorkspace::bxs, nnlsWork_, MahiNnlsWorkspace::nP, MahiNnlsWorkspace::pulseDerivMat, 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));
  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));
  ++nnlsWork_.nP;
}

void MahiFit::onePulseMinimize ( ) const

private

Definition at line 455 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 41 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);
    }
  }

  useTriple=false;
  if(tstrig >= ts4Thresh_ && tsTOT > 0) {
    // 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 518 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 485 of file MahiFit.cc.

References cntsetPulseShape_, and psfPtr_.

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));

}

void MahiFit::resetWorkspace ( ) const

private

Definition at line 660 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::dt, MahiNnlsWorkspace::fullTSOffset, MahiNnlsWorkspace::maxoffset, nnlsWork_, MahiNnlsWorkspace::noiseTerms, MahiNnlsWorkspace::nPulseTot, MahiNnlsWorkspace::pedConstraint, 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;

  nnlsWork_.amplitudes.setZero();
  nnlsWork_.noiseTerms.setZero();
  nnlsWork_.pedConstraint.setZero();

}

void MahiFit::setParameters	(	bool	iDynamicPed,
		double	iTS4Thresh,
		double	chiSqSwitch,
		bool	iApplyTimeSlew,
		HcalTimeSlew::BiasSetting	slewFlavor,
		bool	iCalculateArrivalTime,
		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_, calculateArrivalTime_, 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;

  calculateArrivalTime_ = iCalculateArrivalTime;
  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 472 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 591 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 313 of file MahiFit.cc.

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

Referenced by minimize().

                              {

  SampleMatrix invCovMat;
  invCovMat.setZero(nnlsWork_.tsSize, nnlsWork_.tsSize);
  invCovMat = nnlsWork_.noiseTerms.asDiagonal();
  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 { 
      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(invCovMat);
}

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

private

Definition at line 253 of file MahiFit.cc.

References applyTimeSlew_, HcalTimeSlew::delay(), delta, MahiNnlsWorkspace::dt, hcalTimeSlewDelay_, MahiNnlsWorkspace::maxoffset, meanTime_, nnlsWork_, psfPtr_, slewFlavor_, genVertex_cff::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(float(itQ),slewFlavor_);
  }

  std::array<double, MaxSVSize> pulseN;
  std::array<double, MaxSVSize> pulseM;
  std::array<double, MaxSVSize> 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};

  psfPtr_->singlePulseShapeFuncMahi(&xx[0]);
  psfPtr_->getPulseShape(pulseN);

  psfPtr_->singlePulseShapeFuncMahi(&xxm[0]);
  psfPtr_->getPulseShape(pulseM);
  
  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.5 / nnlsWork_.dt;

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

    pulseShape.coeffRef(iTS+nnlsWork_.maxoffset) = pulseN[iTS+delta];
    pulseDeriv.coeffRef(iTS+nnlsWork_.maxoffset) = (pulseM[iTS+delta]-pulseP[iTS+delta])*invDt;

    pulseM[iTS] -= pulseN[iTS];
    pulseP[iTS] -= pulseN[iTS];
  }

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

Member Data Documentation

std::vector<int> MahiFit::activeBXs_

private

Definition at line 173 of file MahiFit.h.

Referenced by doFit(), and setParameters().

bool MahiFit::applyTimeSlew_

private

Definition at line 164 of file MahiFit.h.

Referenced by setParameters(), and updatePulseShape().

int MahiFit::bxOffsetConf_

private

Definition at line 182 of file MahiFit.h.

Referenced by doFit(), and setParameters().

unsigned int MahiFit::bxSizeConf_

private

Definition at line 181 of file MahiFit.h.

Referenced by doFit(), and setParameters().

bool MahiFit::calculateArrivalTime_

private

Definition at line 168 of file MahiFit.h.

Referenced by doFit(), and setParameters().

float MahiFit::chiSqSwitch_

private

Definition at line 162 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

int MahiFit::cntsetPulseShape_

private

Definition at line 185 of file MahiFit.h.

Referenced by resetPulseShapeTemplate().

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

Definition at line 125 of file MahiFit.h.

Referenced by setPulseShapeTemplate().

float MahiFit::deltaChiSqThresh_

private

Definition at line 178 of file MahiFit.h.

Referenced by minimize(), and setParameters().

bool MahiFit::dynamicPed_

private

Definition at line 160 of file MahiFit.h.

Referenced by doFit(), and setParameters().

const unsigned int MahiFit::fullTSofInterest_

private

Definition at line 151 of file MahiFit.h.

Referenced by phase1Apply().

const unsigned int MahiFit::fullTSSize_

private

Definition at line 150 of file MahiFit.h.

const HcalTimeSlew* MahiFit::hcalTimeSlewDelay_ = 0

Definition at line 126 of file MahiFit.h.

Referenced by setPulseShapeTemplate(), and updatePulseShape().

float MahiFit::meanTime_

private

Definition at line 169 of file MahiFit.h.

Referenced by setParameters(), and updatePulseShape().

int MahiFit::nMaxItersMin_

private

Definition at line 175 of file MahiFit.h.

Referenced by minimize(), and setParameters().

int MahiFit::nMaxItersNNLS_

private

Definition at line 176 of file MahiFit.h.

Referenced by nnls(), and setParameters().

float MahiFit::nnlsThresh_

private

Definition at line 179 of file MahiFit.h.

Referenced by nnls(), and setParameters().

MahiNnlsWorkspace MahiFit::nnlsWork_

mutableprivate

Definition at line 147 of file MahiFit.h.

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

int MahiFit::pedestalBX_ = 100

staticprivate

Definition at line 153 of file MahiFit.h.

Referenced by doFit(), phase1Debug(), and updateCov().

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

private

Definition at line 187 of file MahiFit.h.

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

private

Definition at line 186 of file MahiFit.h.

Referenced by resetPulseShapeTemplate(), and updatePulseShape().

HcalTimeSlew::BiasSetting MahiFit::slewFlavor_

private

Definition at line 165 of file MahiFit.h.

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

float MahiFit::timeLimit_ = 12.5

staticprivate

Definition at line 157 of file MahiFit.h.

Referenced by calculateArrivalTime().

float MahiFit::timeSigmaHPD_

private

Definition at line 170 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

float MahiFit::timeSigmaSiPM_

private

Definition at line 171 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

float MahiFit::ts4Thresh_

private

Definition at line 161 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

float MahiFit::tsDelay1GeV_ =0.f

private

Definition at line 166 of file MahiFit.h.

Referenced by setPulseShapeTemplate(), and updatePulseShape().