Skip to content

Commit

Permalink
Merge pull request #22293 from mariadalfonso/LocalRecoCPUandHLT
Browse files Browse the repository at this point in the history 
HBHE: Local reco cpu
  • Loading branch information
@cmsbuild
cmsbuild committed Mar 2, 2018
2 parents eced59f + c857b1d commit 6ebe5c4
Show file tree
Hide file tree
Showing 7 changed files with 77 additions and 99 deletions.
15 changes: 7 additions & 8 deletions RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h
View file
Expand Up @@ -97,16 +97,16 @@ class MahiFit
float& reconstructedEnergy,
float& reconstructedTime,
bool& useTriple,
float& chi2,
const HcalTimeSlew* hcalTimeSlew_delay) const;
float& chi2) const;

void doFit(std::array<float,3> &correctedOutput, int nbx, const HcalTimeSlew* hcalTimeSlew_delay) const;
void doFit(std::array<float,3> &correctedOutput, const int nbx) const;

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

typedef BXVector::Index Index;
const HcalPulseShapes::Shape* currentPulseShape_=nullptr;
const HcalTimeSlew* hcalTimeSlewDelay_=nullptr;

private:

Expand All @@ -115,8 +115,8 @@ class MahiFit
void updateCov() const;
void updatePulseShape(double itQ, FullSampleVector &pulseShape,
FullSampleVector &pulseDeriv,
FullSampleMatrix &pulseCov,
const HcalTimeSlew* hcalTimeSlew_delay) const;
FullSampleMatrix &pulseCov) const;

double calculateArrivalTime() const;
double calculateChiSq() const;
void nnls() const;
Expand All @@ -127,8 +127,6 @@ class MahiFit

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

double getSiPMDarkCurrent(double darkCurrent, double fcByPE, double lambda) const;

mutable MahiNnlsWorkspace nnlsWork_;

//hard coded in initializer
Expand All @@ -144,6 +142,7 @@ class MahiFit

bool applyTimeSlew_;
HcalTimeSlew::BiasSetting slewFlavor_;
double tsDelay1GeV_=0;

float meanTime_;
float timeSigmaHPD_;
Expand Down
11 changes: 5 additions & 6 deletions RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFitOOTPileupCorrection.h
View file
Expand Up @@ -30,8 +30,7 @@ class PulseShapeFitOOTPileupCorrection
float& reconstructedEnergy,
float& reconstructedTime,
bool & useTriple,
float& chi2,
const HcalTimeSlew* hcalTimeSlew_delay) const;
float& chi2) const;

void setPUParams(bool iPedestalConstraint, bool iTimeConstraint,bool iAddPulseJitter,bool iApplyTimeSlew,
double iTS4Min, const std::vector<double> & iTS4Max,
Expand All @@ -42,16 +41,16 @@ class PulseShapeFitOOTPileupCorrection
const std::vector<double> & its4Chi2, HcalTimeSlew::BiasSetting slewFlavor, int iFitTimes);

const HcalPulseShapes::Shape* currentPulseShape_=nullptr;
const HcalTimeSlew* hcalTimeSlewDelay_=nullptr;
double tsDelay1GeV_=0;

void setPulseShapeTemplate (const HcalPulseShapes::Shape& ps, bool isHPD, unsigned nSamples);
void setPulseShapeTemplate (const HcalPulseShapes::Shape& ps, bool isHPD, unsigned nSamples, const HcalTimeSlew* hcalTimeSlewDelay);
void resetPulseShapeTemplate(const HcalPulseShapes::Shape& ps, unsigned nSamples);

private:

double getSiPMDarkCurrent(double darkCurrent, double fcByPE, double lambda) const;

int pulseShapeFit(const double * energyArr, const double * pedenArr, const double *chargeArr,
const double *pedArr, const double *gainArr, const double tsTOTen, std::vector<float> &fitParsVec, const double * ADCnoise, unsigned int soi, const HcalTimeSlew* hcalTimeSlew_delay) const;
const double *pedArr, const double *gainArr, const double tsTOTen, std::vector<float> &fitParsVec, const double * ADCnoise, unsigned int soi) const;
void fit(int iFit,float &timevalfit,float &chargevalfit,float &pedvalfit,float &chi2,bool &fitStatus,double &iTSMax,
const double &iTSTOTen,double *iEnArr,unsigned (&iBX)[3]) const;

Expand Down
4 changes: 2 additions & 2 deletions RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFunctor.h
View file
Expand Up @@ -23,7 +23,7 @@ namespace FitterFuncs{
unsigned int nSamplesToFit);
~PulseShapeFunctor();

double EvalPulse(const double *pars, unsigned int nPar);
double EvalPulse(const double *pars, const unsigned nPar);

void setDefaultcntNANinfit(){ cntNANinfit =0; }
int getcntNANinfit(){ return cntNANinfit; }
Expand Down Expand Up @@ -52,7 +52,7 @@ namespace FitterFuncs{
std::vector<float> acc25nsVec, diff25nsItvlVec;
std::vector<float> accVarLenIdxZEROVec, diffVarItvlIdxZEROVec;
std::vector<float> accVarLenIdxMinusOneVec, diffVarItvlIdxMinusOneVec;
void funcShape(std::array<double,HcalConst::maxSamples> & ntmpbin, const double &pulseTime, const double &pulseHeight,const double &slew);
void funcShape(std::array<double,HcalConst::maxSamples> & ntmpbin, const double pulseTime, const double pulseHeight,const double slew);
double psFit_x[HcalConst::maxSamples], psFit_y[HcalConst::maxSamples], psFit_erry[HcalConst::maxSamples], psFit_erry2[HcalConst::maxSamples], psFit_slew[HcalConst::maxSamples];

unsigned nSamplesToFit_;
Expand Down
71 changes: 28 additions & 43 deletions RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc
View file
Expand Up @@ -6,11 +6,6 @@ MahiFit::MahiFit() :
fullTSofInterest_(8)
{}

double MahiFit::getSiPMDarkCurrent(double darkCurrent, double fcByPE, double lambda) const {
double mu = darkCurrent * 25 / fcByPE;
return sqrt(mu/pow(1-lambda,3)) * fcByPE;
}

void MahiFit::setParameters(bool iDynamicPed, double iTS4Thresh, double chiSqSwitch,
bool iApplyTimeSlew, HcalTimeSlew::BiasSetting slewFlavor,
double iMeanTime, double iTimeSigmaHPD, double iTimeSigmaSiPM,
Expand Down Expand Up @@ -46,8 +41,7 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,
float& reconstructedEnergy,
float& reconstructedTime,
bool& useTriple,
float& chi2,
const HcalTimeSlew* hcalTimeSlew_delay) const {
float& chi2) const {

assert(channelData.nSamples()==8||channelData.nSamples()==10);

Expand All @@ -61,10 +55,6 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,
if (channelData.hasTimeInfo()) nnlsWork_.dt=timeSigmaSiPM_;
else nnlsWork_.dt=timeSigmaHPD_;

//Dark current value for this channel (SiPM only)
float darkCurrent = getSiPMDarkCurrent(channelData.darkCurrent(),
channelData.fcByPE(),
channelData.lambda());

//Average pedestal width (for covariance matrix constraint)
float pedVal = 0.25*( channelData.tsPedestalWidth(0)*channelData.tsPedestalWidth(0)+
Expand All @@ -88,12 +78,6 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,
//ADC granularity
double noiseADC = (1./sqrt(12))*channelData.tsDFcPerADC(iTS);

//Dark current (for SiPMs)
double noiseDC=0;
if(channelData.hasTimeInfo() && !channelData.hasEffectivePedestals() && (charge-ped)>channelData.tsPedestalWidth(iTS)) {
noiseDC = darkCurrent;
}

//Photostatistics
double noisePhoto = 0;
if ( (charge-ped)>channelData.tsPedestalWidth(iTS)) {
Expand All @@ -104,7 +88,7 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,
double pedWidth = channelData.tsPedestalWidth(iTS);

//Total uncertainty from all sources
nnlsWork_.noiseTerms.coeffRef(iTS) = noiseADC*noiseADC + noiseDC*noiseDC + noisePhoto*noisePhoto + pedWidth*pedWidth;
nnlsWork_.noiseTerms.coeffRef(iTS) = noiseADC*noiseADC + noisePhoto*noisePhoto + pedWidth*pedWidth;

tsTOT += (charge - ped)*channelData.tsGain(0);
if( iTS==nnlsWork_.tsOffset ){
Expand All @@ -118,14 +102,14 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,

// only do pre-fit with 1 pulse if chiSq threshold is positive
if (chiSqSwitch_>0) {
doFit(reconstructedVals,1,hcalTimeSlew_delay);
doFit(reconstructedVals,1);
if (reconstructedVals[2]>chiSqSwitch_) {
doFit(reconstructedVals,0,hcalTimeSlew_delay); //nbx=0 means use configured BXs
doFit(reconstructedVals,0); //nbx=0 means use configured BXs
useTriple=true;
}
}
else {
doFit(reconstructedVals,0,hcalTimeSlew_delay);
doFit(reconstructedVals,0);
useTriple=true;
}
}
Expand All @@ -141,7 +125,7 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,

}

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

unsigned int bxSize=1;

Expand All @@ -159,7 +143,7 @@ void MahiFit::doFit(std::array<float,3> &correctedOutput, int nbx, const HcalTim
nnlsWork_.bxs.coeffRef(0) = 0;
}
else {
for (unsigned int iBX=0; iBX<bxSize; iBX++) {
for (unsigned int iBX=0; iBX<bxSize; ++iBX) {
nnlsWork_.bxs.coeffRef(iBX) = activeBXs_[iBX];
}
}
Expand All @@ -177,7 +161,7 @@ void MahiFit::doFit(std::array<float,3> &correctedOutput, int nbx, const HcalTim
nnlsWork_.errVec = PulseVector::Zero(nnlsWork_.nPulseTot);

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

nnlsWork_.pulseShapeArray[iBX] = FullSampleVector::Zero(MaxFSVSize);
Expand All @@ -192,8 +176,7 @@ void MahiFit::doFit(std::array<float,3> &correctedOutput, int nbx, const HcalTim
updatePulseShape(nnlsWork_.amplitudes.coeff(nnlsWork_.tsOffset + offset),
nnlsWork_.pulseShapeArray[iBX],
nnlsWork_.pulseDerivArray[iBX],
nnlsWork_.pulseCovArray[iBX],
hcalTimeSlew_delay);
nnlsWork_.pulseCovArray[iBX]);

nnlsWork_.ampVec.coeffRef(iBX)=0;

Expand Down Expand Up @@ -233,14 +216,10 @@ void MahiFit::doFit(std::array<float,3> &correctedOutput, int nbx, const HcalTim

double MahiFit::minimize() const {

int iter = 0;
double oldChiSq=9999;
double chiSq=oldChiSq;

while (true) {
if (iter>=nMaxItersMin_) {
break;
}
for( int iter=1; iter<nMaxItersMin_ ; ++iter) {

updateCov();

Expand All @@ -262,20 +241,22 @@ double MahiFit::minimize() const {

if (std::abs(deltaChiSq)<deltaChiSqThresh_) break;

iter++;

}

return chiSq;

}

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

float t0=meanTime_;
if (applyTimeSlew_) t0+=hcalTimeSlew_delay->delay(std::max(1.0, itQ), slewFlavor_);

if(applyTimeSlew_) {
if(itQ<=1.0) t0+=tsDelay1GeV_;
else t0+=hcalTimeSlewDelay_->delay(itQ,slewFlavor_);
}


nnlsWork_.pulseN.fill(0);
nnlsWork_.pulseM.fill(0);
Expand All @@ -298,7 +279,7 @@ void MahiFit::updatePulseShape(double itQ, FullSampleVector &pulseShape, FullSam
//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++) {
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);
Expand All @@ -307,8 +288,8 @@ void MahiFit::updatePulseShape(double itQ, FullSampleVector &pulseShape, FullSam
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++) {
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] );
Expand All @@ -328,7 +309,7 @@ void MahiFit::updateCov() const {
nnlsWork_.invCovMat = nnlsWork_.noiseTerms.asDiagonal();
nnlsWork_.invCovMat +=nnlsWork_.pedConstraint;

for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; iBX++) {
for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
if (nnlsWork_.ampVec.coeff(iBX)==0) continue;

int offset=nnlsWork_.bxs.coeff(iBX);
Expand All @@ -349,7 +330,7 @@ double MahiFit::calculateArrivalTime() const {

int itIndex=0;

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

Expand All @@ -370,7 +351,7 @@ double MahiFit::calculateArrivalTime() const {
void MahiFit::nnls() const {
const unsigned int npulse = nnlsWork_.nPulseTot;

for (unsigned int iBX=0; iBX<npulse; iBX++) {
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);
Expand Down Expand Up @@ -484,10 +465,14 @@ double MahiFit::calculateChiSq() const {
return (nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat*nnlsWork_.ampVec - nnlsWork_.amplitudes)).squaredNorm();
}

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

if (!(&ps == currentPulseShape_ ))
{

hcalTimeSlewDelay_ = hcalTimeSlewDelay;
tsDelay1GeV_= hcalTimeSlewDelay->delay(1.0, slewFlavor_);

resetPulseShapeTemplate(ps);
currentPulseShape_ = &ps;
}
Expand Down

0 comments on commit 6ebe5c4

Please sign in to comment.