CMS 3D CMS Logo

MahiFit.cc
Go to the documentation of this file.
3 
5  fullTSSize_(19),
6  fullTSofInterest_(8)
7 {}
8 
9 void MahiFit::setParameters(bool iDynamicPed, double iTS4Thresh, double chiSqSwitch,
10  bool iApplyTimeSlew, HcalTimeSlew::BiasSetting slewFlavor,
11  double iMeanTime, double iTimeSigmaHPD, double iTimeSigmaSiPM,
12  const std::vector <int> &iActiveBXs, int iNMaxItersMin, int iNMaxItersNNLS,
13  double iDeltaChiSqThresh, double iNnlsThresh) {
14 
15  dynamicPed_ = iDynamicPed;
16 
17  ts4Thresh_ = iTS4Thresh;
19 
20  applyTimeSlew_ = iApplyTimeSlew;
21  slewFlavor_ = slewFlavor;
22 
23  meanTime_ = iMeanTime;
24  timeSigmaHPD_ = iTimeSigmaHPD;
25  timeSigmaSiPM_ = iTimeSigmaSiPM;
26 
27  activeBXs_ = iActiveBXs;
28 
29  nMaxItersMin_ = iNMaxItersMin;
30  nMaxItersNNLS_ = iNMaxItersNNLS;
31 
32  deltaChiSqThresh_ = iDeltaChiSqThresh;
33  nnlsThresh_ = iNnlsThresh;
34 
35  bxOffsetConf_ = -(*std::min_element(activeBXs_.begin(), activeBXs_.end()));
36  bxSizeConf_ = activeBXs_.size();
37 
38 }
39 
40 void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,
41  float& reconstructedEnergy,
42  float& reconstructedTime,
43  bool& useTriple,
44  float& chi2) const {
45 
46  assert(channelData.nSamples()==8||channelData.nSamples()==10);
47 
49 
50  nnlsWork_.tsSize = channelData.nSamples();
51  nnlsWork_.tsOffset = channelData.soi();
53 
54  // 1 sigma time constraint
55  if (channelData.hasTimeInfo()) nnlsWork_.dt=timeSigmaSiPM_;
57 
58 
59  //Average pedestal width (for covariance matrix constraint)
60  float pedVal = 0.25*( channelData.tsPedestalWidth(0)*channelData.tsPedestalWidth(0)+
61  channelData.tsPedestalWidth(1)*channelData.tsPedestalWidth(1)+
62  channelData.tsPedestalWidth(2)*channelData.tsPedestalWidth(2)+
63  channelData.tsPedestalWidth(3)*channelData.tsPedestalWidth(3) );
64 
65  nnlsWork_.pedConstraint = pedVal*SampleMatrix::Ones(nnlsWork_.tsSize, nnlsWork_.tsSize);
68 
69  std::array<float,3> reconstructedVals {{ 0.0, -9999, -9999 }};
70 
71  double tsTOT = 0, tstrig = 0; // in GeV
72  for(unsigned int iTS=0; iTS<nnlsWork_.tsSize; ++iTS){
73  double charge = channelData.tsRawCharge(iTS);
74  double ped = channelData.tsPedestal(iTS);
75 
76  nnlsWork_.amplitudes.coeffRef(iTS) = charge - ped;
77 
78  //ADC granularity
79  double noiseADC = (1./sqrt(12))*channelData.tsDFcPerADC(iTS);
80 
81  //Photostatistics
82  double noisePhoto = 0;
83  if ( (charge-ped)>channelData.tsPedestalWidth(iTS)) {
84  noisePhoto = sqrt((charge-ped)*channelData.fcByPE());
85  }
86 
87  //Electronic pedestal
88  double pedWidth = channelData.tsPedestalWidth(iTS);
89 
90  //Total uncertainty from all sources
91  nnlsWork_.noiseTerms.coeffRef(iTS) = noiseADC*noiseADC + noisePhoto*noisePhoto + pedWidth*pedWidth;
92 
93  tsTOT += (charge - ped)*channelData.tsGain(0);
94  if( iTS==nnlsWork_.tsOffset ){
95  tstrig += (charge - ped)*channelData.tsGain(0);
96  }
97  }
98 
99  if(tstrig >= ts4Thresh_ && tsTOT > 0) {
100 
101  useTriple=false;
102 
103  // only do pre-fit with 1 pulse if chiSq threshold is positive
104  if (chiSqSwitch_>0) {
105  doFit(reconstructedVals,1);
106  if (reconstructedVals[2]>chiSqSwitch_) {
107  doFit(reconstructedVals,0); //nbx=0 means use configured BXs
108  useTriple=true;
109  }
110  }
111  else {
112  doFit(reconstructedVals,0);
113  useTriple=true;
114  }
115  }
116  else{
117  reconstructedVals.at(0) = 0.; //energy
118  reconstructedVals.at(1) = -9999.; //time
119  reconstructedVals.at(2) = -9999.; //chi2
120  }
121 
122  reconstructedEnergy = reconstructedVals[0]*channelData.tsGain(0);
123  reconstructedTime = reconstructedVals[1];
124  chi2 = reconstructedVals[2];
125 
126 }
127 
128 void MahiFit::doFit(std::array<float,3> &correctedOutput, int nbx) const {
129 
130  unsigned int bxSize=1;
131 
132  if (nbx==1) {
133  nnlsWork_.bxOffset = 0;
134  }
135  else {
136  bxSize = bxSizeConf_;
138  }
139 
140  nnlsWork_.bxs.resize(bxSize);
141 
142  if (nbx==1) {
143  nnlsWork_.bxs.coeffRef(0) = 0;
144  }
145  else {
146  for (unsigned int iBX=0; iBX<bxSize; ++iBX) {
147  nnlsWork_.bxs.coeffRef(iBX) = activeBXs_[iBX];
148  }
149  }
150 
151  nnlsWork_.nPulseTot = bxSize;
152 
153  if (dynamicPed_) {
157  }
158 
160  nnlsWork_.ampVec = PulseVector::Zero(nnlsWork_.nPulseTot);
161  nnlsWork_.errVec = PulseVector::Zero(nnlsWork_.nPulseTot);
162 
163  int offset=0;
164  for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
165  offset=nnlsWork_.bxs.coeff(iBX);
166 
167  nnlsWork_.pulseShapeArray[iBX] = FullSampleVector::Zero(MaxFSVSize);
168  nnlsWork_.pulseDerivArray[iBX] = FullSampleVector::Zero(MaxFSVSize);
170 
171  if (offset==pedestalBX_) {
172  nnlsWork_.ampVec.coeffRef(iBX) = 0;
173  }
174  else {
175 
179  nnlsWork_.pulseCovArray[iBX]);
180 
181  nnlsWork_.ampVec.coeffRef(iBX)=0;
182 
184 
185  }
186  }
187 
188  nnlsWork_.pulseMat.col(nnlsWork_.nPulseTot-1) = SampleVector::Ones(nnlsWork_.tsSize);
189 
192 
193  double chiSq = minimize();
194 
196 
197  bool foundintime = false;
198  unsigned int ipulseintime = 0;
199 
200  for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
201  if (nnlsWork_.bxs.coeff(iBX)==0) {
202  ipulseintime = iBX;
203  foundintime = true;
204  }
205  }
206 
207  if (foundintime) {
208  correctedOutput.at(0) = nnlsWork_.ampVec.coeff(ipulseintime); //charge
209  double arrivalTime = calculateArrivalTime();
210  correctedOutput.at(1) = arrivalTime; //time
211  correctedOutput.at(2) = chiSq; //chi2
212 
213  }
214 
215 }
216 
217 double MahiFit::minimize() const {
218 
219  double oldChiSq=9999;
220  double chiSq=oldChiSq;
221 
222  for( int iter=1; iter<nMaxItersMin_ ; ++iter) {
223 
224  updateCov();
225 
226  if (nnlsWork_.nPulseTot>1) {
227  nnls();
228  }
229  else {
231  }
232 
233  double newChiSq=calculateChiSq();
234  double deltaChiSq = newChiSq - chiSq;
235 
236  if (newChiSq==oldChiSq && newChiSq<chiSq) {
237  break;
238  }
239  oldChiSq=chiSq;
240  chiSq = newChiSq;
241 
242  if (std::abs(deltaChiSq)<deltaChiSqThresh_) break;
243 
244  }
245 
246  return chiSq;
247 
248 }
249 
250 void MahiFit::updatePulseShape(double itQ, FullSampleVector &pulseShape, FullSampleVector &pulseDeriv,
251  FullSampleMatrix &pulseCov) const {
252 
253  float t0=meanTime_;
254 
255  if(applyTimeSlew_) {
256  if(itQ<=1.0) t0+=tsDelay1GeV_;
257  else t0+=hcalTimeSlewDelay_->delay(itQ,slewFlavor_);
258  }
259 
260 
261  nnlsWork_.pulseN.fill(0);
262  nnlsWork_.pulseM.fill(0);
263  nnlsWork_.pulseP.fill(0);
264 
265  const double xx[4]={t0, 1.0, 0.0, 3};
266  const double xxm[4]={-nnlsWork_.dt+t0, 1.0, 0.0, 3};
267  const double xxp[4]={ nnlsWork_.dt+t0, 1.0, 0.0, 3};
268 
269  (*pfunctor_)(&xx[0]);
270  psfPtr_->getPulseShape(nnlsWork_.pulseN);
271 
272  (*pfunctor_)(&xxm[0]);
273  psfPtr_->getPulseShape(nnlsWork_.pulseM);
274 
275  (*pfunctor_)(&xxp[0]);
276  psfPtr_->getPulseShape(nnlsWork_.pulseP);
277 
278  //in the 2018+ case where the sample of interest (SOI) is in TS3, add an extra offset to align
279  //with previous SOI=TS4 case assumed by psfPtr_->getPulseShape()
280  int delta =nnlsWork_. tsOffset == 3 ? 1 : 0;
281 
282  for (unsigned int iTS=nnlsWork_.fullTSOffset; iTS<nnlsWork_.fullTSOffset + nnlsWork_.tsSize; ++iTS) {
283 
284  pulseShape.coeffRef(iTS) = nnlsWork_.pulseN[iTS-nnlsWork_.fullTSOffset+delta];
285  pulseDeriv.coeffRef(iTS) = 0.5*(nnlsWork_.pulseM[iTS-nnlsWork_.fullTSOffset+delta]+nnlsWork_.pulseP[iTS-nnlsWork_.fullTSOffset+delta])/(2*nnlsWork_.dt);
286 
289  }
290 
291  for (unsigned int iTS=nnlsWork_.fullTSOffset; iTS<nnlsWork_.fullTSOffset+nnlsWork_.tsSize; ++iTS) {
292  for (unsigned int jTS=nnlsWork_.fullTSOffset; jTS<iTS+1; ++jTS) {
293 
296 
297  pulseCov(iTS,jTS) += tmp;
298  pulseCov(jTS,iTS) += tmp;
299 
300  }
301  }
302 
303 }
304 
305 void MahiFit::updateCov() const {
306 
308 
309  nnlsWork_.invCovMat = nnlsWork_.noiseTerms.asDiagonal();
311 
312  for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
313  if (nnlsWork_.ampVec.coeff(iBX)==0) continue;
314 
315  int offset=nnlsWork_.bxs.coeff(iBX);
316 
317  if (offset==pedestalBX_) continue;
318  else {
319  nnlsWork_.invCovMat += nnlsWork_.ampVec.coeff(iBX)*nnlsWork_.ampVec.coeff(iBX)
321  }
322  }
323 
325 }
326 
328 
330 
331  int itIndex=0;
332 
333  for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
334  int offset=nnlsWork_.bxs.coeff(iBX);
335  if (offset==0) itIndex=iBX;
336 
337  if (offset==pedestalBX_) {
338  nnlsWork_.pulseDerivMat.col(iBX) = SampleVector::Zero(nnlsWork_.tsSize);
339  }
340  else {
342  }
343  }
344 
345  PulseVector solution = nnlsWork_.pulseDerivMat.colPivHouseholderQr().solve(nnlsWork_.residuals);
346  return solution.coeff(itIndex)/nnlsWork_.ampVec.coeff(itIndex);
347 
348 }
349 
350 
351 void MahiFit::nnls() const {
352  const unsigned int npulse = nnlsWork_.nPulseTot;
353 
354  for (unsigned int iBX=0; iBX<npulse; ++iBX) {
355  int offset=nnlsWork_.bxs.coeff(iBX);
356  if (offset==pedestalBX_) {
357  nnlsWork_.pulseMat.col(iBX) = SampleVector::Ones(nnlsWork_.tsSize);
358  }
359  else {
361  }
362  }
363 
365  nnlsWork_.aTaMat = nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.invcovp);
366  nnlsWork_.aTbVec = nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes));
367 
368  int iter = 0;
369  Index idxwmax = 0;
370  double wmax = 0.0;
371  double threshold = nnlsThresh_;
372 
373  nnlsWork_.nP=0;
374 
375  while (true) {
376  if (iter>0 || nnlsWork_.nP==0) {
377  if ( nnlsWork_.nP==std::min(npulse, nnlsWork_.tsSize)) break;
378 
379  const unsigned int nActive = npulse - nnlsWork_.nP;
381 
382  Index idxwmaxprev = idxwmax;
383  double wmaxprev = wmax;
384  wmax = nnlsWork_.updateWork.tail(nActive).maxCoeff(&idxwmax);
385 
386  if (wmax<threshold || (idxwmax==idxwmaxprev && wmax==wmaxprev)) {
387  break;
388  }
389 
390  if (iter>=nMaxItersNNLS_) {
391  break;
392  }
393 
394  //unconstrain parameter
395  Index idxp = nnlsWork_.nP + idxwmax;
397 
398  }
399 
400  while (true) {
401  if (nnlsWork_.nP==0) break;
402 
404 
406 
407  //check solution
408  bool positive = true;
409  for (unsigned int i = 0; i < nnlsWork_.nP; ++i)
410  positive &= (nnlsWork_.ampvecpermtest(i) > 0);
411  if (positive) {
413  break;
414  }
415 
416  //update parameter vector
417  Index minratioidx=0;
418 
419  // no realizable optimization here (because it autovectorizes!)
420  double minratio = std::numeric_limits<double>::max();
421  for (unsigned int ipulse=0; ipulse<nnlsWork_.nP; ++ipulse) {
422  if (nnlsWork_.ampvecpermtest.coeff(ipulse)<=0.) {
423  const double c_ampvec = nnlsWork_.ampVec.coeff(ipulse);
424  const double ratio = c_ampvec/(c_ampvec-nnlsWork_.ampvecpermtest.coeff(ipulse));
425  if (ratio<minratio) {
426  minratio = ratio;
427  minratioidx = ipulse;
428  }
429  }
430  }
432 
433  //avoid numerical problems with later ==0. check
434  nnlsWork_.ampVec.coeffRef(minratioidx) = 0.;
435 
436  nnlsConstrainParameter(minratioidx);
437  }
438 
439  ++iter;
440 
441  //adaptive convergence threshold to avoid infinite loops but still
442  //ensure best value is used
443  if (iter%10==0) {
444  threshold *= 10.;
445  }
446 
447  }
448 
449 }
450 
452 
454 
455  SingleMatrix aTamatval = nnlsWork_.invcovp.transpose()*nnlsWork_.invcovp;
456  SingleVector aTbvecval = nnlsWork_.invcovp.transpose()*nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes);
457 
458  nnlsWork_.ampVec.coeffRef(0) = std::max(0., aTbvecval.coeff(0)/aTamatval.coeff(0));
459 
460 
461 }
462 
463 double MahiFit::calculateChiSq() const {
464 
465  return (nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat*nnlsWork_.ampVec - nnlsWork_.amplitudes)).squaredNorm();
466 }
467 
468 void MahiFit::setPulseShapeTemplate(const HcalPulseShapes::Shape& ps,const HcalTimeSlew* hcalTimeSlewDelay) {
469 
470  if (!(&ps == currentPulseShape_ ))
471  {
472 
473  hcalTimeSlewDelay_ = hcalTimeSlewDelay;
474  tsDelay1GeV_= hcalTimeSlewDelay->delay(1.0, slewFlavor_);
475 
477  currentPulseShape_ = &ps;
478  }
479 }
480 
483 
484  // only the pulse shape itself from PulseShapeFunctor is used for Mahi
485  // the uncertainty terms calculated inside PulseShapeFunctor are used for Method 2 only
486  psfPtr_.reset(new FitterFuncs::PulseShapeFunctor(ps,false,false,false,
487  1,0,0,10));
488  pfunctor_ = std::unique_ptr<ROOT::Math::Functor>( new ROOT::Math::Functor(psfPtr_.get(),&FitterFuncs::PulseShapeFunctor::singlePulseShapeFunc, 3) );
489 
490 
491 }
492 
494  nnlsWork_.aTaMat.col(nnlsWork_.nP).swap(nnlsWork_.aTaMat.col(idxp));
495  nnlsWork_.aTaMat.row(nnlsWork_.nP).swap(nnlsWork_.aTaMat.row(idxp));
496  nnlsWork_.pulseMat.col(nnlsWork_.nP).swap(nnlsWork_.pulseMat.col(idxp));
497  std::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP),nnlsWork_.aTbVec.coeffRef(idxp));
498  std::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP),nnlsWork_.ampVec.coeffRef(idxp));
499  std::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP),nnlsWork_.bxs.coeffRef(idxp));
500  ++nnlsWork_.nP;
501 }
502 
503 void MahiFit::nnlsConstrainParameter(Index minratioidx) const {
504  nnlsWork_.aTaMat.col(nnlsWork_.nP-1).swap(nnlsWork_.aTaMat.col(minratioidx));
505  nnlsWork_.aTaMat.row(nnlsWork_.nP-1).swap(nnlsWork_.aTaMat.row(minratioidx));
506  nnlsWork_.pulseMat.col(nnlsWork_.nP-1).swap(nnlsWork_.pulseMat.col(minratioidx));
507  std::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP-1),nnlsWork_.aTbVec.coeffRef(minratioidx));
508  std::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP-1),nnlsWork_.ampVec.coeffRef(minratioidx));
509  std::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP-1),nnlsWork_.bxs.coeffRef(minratioidx));
510  --nnlsWork_.nP;
511 
512 }
513 
514 void MahiFit::solveSubmatrix(PulseMatrix& mat, PulseVector& invec, PulseVector& outvec, unsigned nP) const {
515  using namespace Eigen;
516  switch( nP ) { // pulse matrix is always square.
517  case 10:
518  {
519  Matrix<double,10,10> temp = mat;
520  outvec.head<10>() = temp.ldlt().solve(invec.head<10>());
521  }
522  break;
523  case 9:
524  {
525  Matrix<double,9,9> temp = mat.topLeftCorner<9,9>();
526  outvec.head<9>() = temp.ldlt().solve(invec.head<9>());
527  }
528  break;
529  case 8:
530  {
531  Matrix<double,8,8> temp = mat.topLeftCorner<8,8>();
532  outvec.head<8>() = temp.ldlt().solve(invec.head<8>());
533  }
534  break;
535  case 7:
536  {
537  Matrix<double,7,7> temp = mat.topLeftCorner<7,7>();
538  outvec.head<7>() = temp.ldlt().solve(invec.head<7>());
539  }
540  break;
541  case 6:
542  {
543  Matrix<double,6,6> temp = mat.topLeftCorner<6,6>();
544  outvec.head<6>() = temp.ldlt().solve(invec.head<6>());
545  }
546  break;
547  case 5:
548  {
549  Matrix<double,5,5> temp = mat.topLeftCorner<5,5>();
550  outvec.head<5>() = temp.ldlt().solve(invec.head<5>());
551  }
552  break;
553  case 4:
554  {
555  Matrix<double,4,4> temp = mat.topLeftCorner<4,4>();
556  outvec.head<4>() = temp.ldlt().solve(invec.head<4>());
557  }
558  break;
559  case 3:
560  {
561  Matrix<double,3,3> temp = mat.topLeftCorner<3,3>();
562  outvec.head<3>() = temp.ldlt().solve(invec.head<3>());
563  }
564  break;
565  case 2:
566  {
567  Matrix<double,2,2> temp = mat.topLeftCorner<2,2>();
568  outvec.head<2>() = temp.ldlt().solve(invec.head<2>());
569  }
570  break;
571  case 1:
572  {
573  Matrix<double,1,1> temp = mat.topLeftCorner<1,1>();
574  outvec.head<1>() = temp.ldlt().solve(invec.head<1>());
575  }
576  break;
577  default:
578  throw cms::Exception("HcalMahiWeirdState")
579  << "Weird number of pulses encountered in Mahi, module is configured incorrectly!";
580  }
581 }
582 
584 
586  nnlsWork_.tsSize=0;
590  nnlsWork_.dt=0;
591 
595 
599 
600  nnlsWork_.amplitudes.setZero();
601  nnlsWork_.bxs.setZero();
602  nnlsWork_.invCovMat.setZero();
603  nnlsWork_.noiseTerms.setZero();
604  nnlsWork_.pedConstraint.setZero();
605  nnlsWork_.pulseMat.setZero();
606  nnlsWork_.pulseDerivMat.setZero();
607  nnlsWork_.residuals.setZero();
608  nnlsWork_.ampVec.setZero();
609  nnlsWork_.errVec.setZero();
610  nnlsWork_.ampvecpermtest.setZero();
611  nnlsWork_.invcovp.setZero();
612  nnlsWork_.aTaMat.setZero();
613  nnlsWork_.aTbVec.setZero();
614  nnlsWork_.updateWork.setZero();
615  nnlsWork_.covDecompLinv.setZero();
616  nnlsWork_.topleft_work.setZero();
617 
618 
619 
620 }
unsigned int nPulseTot
Definition: MahiFit.h:19
dbl * delta
Definition: mlp_gen.cc:36
void nnls() const
Definition: MahiFit.cc:351
void phase1Apply(const HBHEChannelInfo &channelData, float &reconstructedEnergy, float &reconstructedTime, bool &useTriple, float &chi2) const
Definition: MahiFit.cc:40
float chiSqSwitch_
Definition: MahiFit.h:141
SamplePulseMatrix invcovp
Definition: MahiFit.h:72
HcalTimeSlew::BiasSetting slewFlavor_
Definition: MahiFit.h:144
bool hasTimeInfo() const
double tsGain(const unsigned ts) const
void nnlsConstrainParameter(Index minratioidx) const
Definition: MahiFit.cc:503
SampleVector amplitudes
Definition: MahiFit.h:30
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:130
BXVector::Index Index
Definition: MahiFit.h:107
double singlePulseShapeFunc(const double *x)
void solveSubmatrix(PulseMatrix &mat, PulseVector &invec, PulseVector &outvec, unsigned nP) const
Definition: MahiFit.cc:514
double delay(double fC, BiasSetting bias=Medium) const
Returns the amount (ns) by which a pulse of the given number of fC will be delayed by the timeslew ef...
Definition: HcalTimeSlew.cc:14
static constexpr int pedestalBX_
Definition: MahiFit.h:136
bool applyTimeSlew_
Definition: MahiFit.h:143
void resetPulseShapeTemplate(const HcalPulseShapes::Shape &ps)
Definition: MahiFit.cc:481
unsigned soi() const
PulseVector ampVec
Definition: MahiFit.h:67
SampleDecompLLT covDecomp
Definition: MahiFit.h:77
PulseVector ampvecpermtest
Definition: MahiFit.h:70
constexpr int MaxFSVSize
std::array< FullSampleVector, MaxPVSize > pulseShapeArray
Definition: MahiFit.h:46
float meanTime_
Definition: MahiFit.h:147
void nnlsUnconstrainParameter(Index idxp) const
Definition: MahiFit.cc:493
std::array< double, MaxSVSize > pulseM
Definition: MahiFit.h:53
PulseMatrix topleft_work
Definition: MahiFit.h:79
double tsDelay1GeV_
Definition: MahiFit.h:145
double tsPedestal(const unsigned ts) const
SamplePulseMatrix pulseDerivMat
Definition: MahiFit.h:60
std::array< double, MaxSVSize > pulseN
Definition: MahiFit.h:52
int bxOffsetConf_
Definition: MahiFit.h:160
PulseVector aTbVec
Definition: MahiFit.h:74
double tsRawCharge(const unsigned ts) const
float nnlsThresh_
Definition: MahiFit.h:157
int cntsetPulseShape_
Definition: MahiFit.h:163
void onePulseMinimize() const
Definition: MahiFit.cc:451
std::unique_ptr< FitterFuncs::PulseShapeFunctor > psfPtr_
Definition: MahiFit.h:164
MahiFit()
Definition: MahiFit.cc:4
Eigen::Matrix< double, FullSampleVectorSize, 1 > FullSampleVector
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
T sqrt(T t)
Definition: SSEVec.h:18
void doFit(std::array< float, 3 > &correctedOutput, const int nbx) const
Definition: MahiFit.cc:128
unsigned int bxSizeConf_
Definition: MahiFit.h:159
double fcByPE() const
int nMaxItersMin_
Definition: MahiFit.h:153
PulseVector residuals
Definition: MahiFit.h:63
Eigen::Matrix< double, 1, 1 > SingleVector
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
#define end
Definition: vmac.h:39
T min(T a, T b)
Definition: MathUtil.h:58
double calculateArrivalTime() const
Definition: MahiFit.cc:327
SampleMatrix pedConstraint
Definition: MahiFit.h:39
std::vector< int > activeBXs_
Definition: MahiFit.h:151
std::array< FullSampleVector, MaxPVSize > pulseDerivArray
Definition: MahiFit.h:49
const unsigned int fullTSofInterest_
Definition: MahiFit.h:134
Polynomial< 0 > Constant
Definition: Constant.h:6
SampleMatrix covDecompLinv
Definition: MahiFit.h:78
SampleMatrix invCovMat
Definition: MahiFit.h:33
void updatePulseShape(double itQ, FullSampleVector &pulseShape, FullSampleVector &pulseDeriv, FullSampleMatrix &pulseCov) const
Definition: MahiFit.cc:250
unsigned int tsOffset
Definition: MahiFit.h:21
void resetWorkspace() const
Definition: MahiFit.cc:583
std::unique_ptr< ROOT::Math::Functor > pfunctor_
Definition: MahiFit.h:165
Eigen::Matrix< double, FullSampleVectorSize, FullSampleVectorSize > FullSampleMatrix
void resize(int bx, unsigned size)
double tsPedestalWidth(const unsigned ts) const
double calculateChiSq() const
Definition: MahiFit.cc:463
PulseMatrix aTaMat
Definition: MahiFit.h:73
unsigned int fullTSOffset
Definition: MahiFit.h:22
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, PulseVectorSize, 1 > PulseVector
SamplePulseMatrix pulseMat
Definition: MahiFit.h:57
bool dynamicPed_
Definition: MahiFit.h:139
std::array< FullSampleMatrix, MaxPVSize > pulseCovArray
Definition: MahiFit.h:43
float timeSigmaSiPM_
Definition: MahiFit.h:149
SampleVector noiseTerms
Definition: MahiFit.h:36
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
#define begin
Definition: vmac.h:32
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)
Definition: MahiFit.cc:9
float ts4Thresh_
Definition: MahiFit.h:140
const HcalTimeSlew * hcalTimeSlewDelay_
Definition: MahiFit.h:109
unsigned int tsSize
Definition: MahiFit.h:20
float deltaChiSqThresh_
Definition: MahiFit.h:156
void updateCov() const
Definition: MahiFit.cc:305
float timeSigmaHPD_
Definition: MahiFit.h:148
unsigned int nP
Definition: MahiFit.h:66
BXVector bxs
Definition: MahiFit.h:27
const HcalPulseShapes::Shape * currentPulseShape_
Definition: MahiFit.h:108
unsigned nSamples() const
PulseVector errVec
Definition: MahiFit.h:69
int nMaxItersNNLS_
Definition: MahiFit.h:154
std::array< double, MaxSVSize > pulseP
Definition: MahiFit.h:54
Eigen::Matrix< double, 1, 1 > SingleMatrix
PulseVector updateWork
Definition: MahiFit.h:75
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, PulseVectorSize, PulseVectorSize > PulseMatrix
double minimize() const
Definition: MahiFit.cc:217
float tsDFcPerADC(const unsigned ts) const
void setPulseShapeTemplate(const HcalPulseShapes::Shape &ps, const HcalTimeSlew *hcalTimeSlewDelay)
Definition: MahiFit.cc:468