EcalRecAlgos/src/EcalUncalibRecHitMultiFitAlgo.cc

0001 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitMultiFitAlgo.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "CondFormats/EcalObjects/interface/EcalPedestals.h"
0006 #include "CondFormats/EcalObjects/interface/EcalGainRatios.h"
0007
0008 EcalUncalibRecHitMultiFitAlgo::EcalUncalibRecHitMultiFitAlgo()
0009     : _computeErrors(true),
0010       _doPrefit(false),
0011       _prefitMaxChiSq(1.0),
0012       _dynamicPedestals(false),
0013       _mitigateBadSamples(false),
0014       _selectiveBadSampleCriteria(false),
0015       _addPedestalUncertainty(0.),
0016       _simplifiedNoiseModelForGainSwitch(true),
0017       _gainSwitchUseMaxSample(false) {
0018   _singlebx.resize(1);
0019   _singlebx << 0;
0020
0021   _pulsefuncSingle.disableErrorCalculation();
0022   _pulsefuncSingle.setMaxIters(1);
0023   _pulsefuncSingle.setMaxIterWarnings(false);
0024 }
0025
0026 /// compute rechits
0027 EcalUncalibratedRecHit EcalUncalibRecHitMultiFitAlgo::makeRecHit(const EcalDataFrame &dataFrame,
0028                                                                  const EcalPedestals::Item *aped,
0029                                                                  const EcalMGPAGainRatio *aGain,
0030                                                                  const SampleMatrixGainArray &noisecors,
0031                                                                  const FullSampleVector &fullpulse,
0032                                                                  const FullSampleMatrix &fullpulsecov,
0033                                                                  const BXVector &activeBX) {
0034   uint32_t flags = 0;
0035
0036   const unsigned int nsample = EcalDataFrame::MAXSAMPLES;
0037
0038   double maxamplitude = -std::numeric_limits<double>::max();
0039   const unsigned int iSampleMax = 5;
0040   const unsigned int iFullPulseMax = 9;
0041
0042   double pedval = 0.;
0043
0044   SampleVector amplitudes;
0045   SampleGainVector gainsNoise;
0046   SampleGainVector gainsPedestal;
0047   SampleGainVector badSamples = SampleGainVector::Zero();
0048   bool hasSaturation = dataFrame.isSaturated();
0049   bool hasGainSwitch = hasSaturation || dataFrame.hasSwitchToGain6() || dataFrame.hasSwitchToGain1();
0050
0051   //no dynamic pedestal in case of gain switch, since then the fit becomes too underconstrained
0052   bool dynamicPedestal = _dynamicPedestals && !hasGainSwitch;
0053
0054   for (unsigned int iSample = 0; iSample < nsample; iSample++) {
0055     const EcalMGPASample &sample = dataFrame.sample(iSample);
0056
0057     double amplitude = 0.;
0058     int gainId = sample.gainId();
0059
0060     double pedestal = 0.;
0061     double gainratio = 1.;
0062
0063     if (gainId == 0 || gainId == 3) {
0064       pedestal = aped->mean_x1;
0065       gainratio = aGain->gain6Over1() * aGain->gain12Over6();
0066       gainsNoise[iSample] = 2;
0067       gainsPedestal[iSample] = dynamicPedestal ? 2 : -1;  //-1 for static pedestal
0068     } else if (gainId == 1) {
0069       pedestal = aped->mean_x12;
0070       gainratio = 1.;
0071       gainsNoise[iSample] = 0;
0072       gainsPedestal[iSample] = dynamicPedestal ? 0 : -1;  //-1 for static pedestal
0073     } else if (gainId == 2) {
0074       pedestal = aped->mean_x6;
0075       gainratio = aGain->gain12Over6();
0076       gainsNoise[iSample] = 1;
0077       gainsPedestal[iSample] = dynamicPedestal ? 1 : -1;  //-1 for static pedestals
0078     }
0079
0080     if (dynamicPedestal) {
0081       amplitude = (double)(sample.adc()) * gainratio;
0082     } else {
0083       amplitude = ((double)(sample.adc()) - pedestal) * gainratio;
0084     }
0085
0086     if (gainId == 0) {
0087       edm::LogError("EcalUncalibRecHitMultiFitAlgo")
0088           << "Saturation encountered.  Multifit is not intended to be used for saturated channels.";
0089       //saturation
0090       if (dynamicPedestal) {
0091         amplitude = 4095. * gainratio;
0092       } else {
0093         amplitude = (4095. - pedestal) * gainratio;
0094       }
0095     }
0096
0097     amplitudes[iSample] = amplitude;
0098
0099     if (iSample == iSampleMax) {
0100       maxamplitude = amplitude;
0101       pedval = pedestal;
0102     }
0103   }
0104
0105   double amplitude, amperr, chisq;
0106   bool status = false;
0107
0108   //special handling for gain switch, where sample before maximum is potentially affected by slew rate limitation
0109   //optionally apply a stricter criteria, assuming slew rate limit is only reached in case where maximum sample has gain switched but previous sample has not
0110   //option 1: use simple max-sample algorithm
0111   if (hasGainSwitch && _gainSwitchUseMaxSample) {
0112     double maxpulseamplitude = maxamplitude / fullpulse[iFullPulseMax];
0113     EcalUncalibratedRecHit rh(dataFrame.id(), maxpulseamplitude, pedval, 0., 0., flags);
0114     rh.setAmplitudeError(0.);
0115     for (unsigned int ipulse = 0; ipulse < _pulsefunc.BXs().rows(); ++ipulse) {
0116       int bx = _pulsefunc.BXs().coeff(ipulse);
0117       if (bx != 0) {
0118         rh.setOutOfTimeAmplitude(bx + 5, 0.0);
0119       }
0120     }
0121     return rh;
0122   }
0123
0124   //option2: A floating negative single-sample offset is added to the fit
0125   //such that the affected sample is treated only as a lower limit for the true amplitude
0126   bool mitigateBadSample = _mitigateBadSamples && hasGainSwitch && iSampleMax > 0;
0127   mitigateBadSample &=
0128       (!_selectiveBadSampleCriteria || (gainsNoise.coeff(iSampleMax - 1) != gainsNoise.coeff(iSampleMax)));
0129   if (mitigateBadSample) {
0130     badSamples[iSampleMax - 1] = 1;
0131   }
0132
0133   //compute noise covariance matrix, which depends on the sample gains
0134   SampleMatrix noisecov;
0135   if (hasGainSwitch) {
0136     std::array<double, 3> pedrmss = {{aped->rms_x12, aped->rms_x6, aped->rms_x1}};
0137     std::array<double, 3> gainratios = {{1., aGain->gain12Over6(), aGain->gain6Over1() * aGain->gain12Over6()}};
0138     if (_simplifiedNoiseModelForGainSwitch) {
0139       int gainidxmax = gainsNoise[iSampleMax];
0140       noisecov = gainratios[gainidxmax] * gainratios[gainidxmax] * pedrmss[gainidxmax] * pedrmss[gainidxmax] *
0141                  noisecors[gainidxmax];
0142       if (!dynamicPedestal && _addPedestalUncertainty > 0.) {
0143         //add fully correlated component to noise covariance to inflate pedestal uncertainty
0144         noisecov += _addPedestalUncertainty * _addPedestalUncertainty * SampleMatrix::Ones();
0145       }
0146     } else {
0147       noisecov = SampleMatrix::Zero();
0148       for (unsigned int gainidx = 0; gainidx < noisecors.size(); ++gainidx) {
0149         SampleGainVector mask = gainidx * SampleGainVector::Ones();
0150         SampleVector pedestal = (gainsNoise.array() == mask.array()).cast<SampleVector::value_type>();
0151         if (pedestal.maxCoeff() > 0.) {
0152           //select out relevant components of each correlation matrix, and assume no correlation between samples with
0153           //different gain
0154           noisecov += gainratios[gainidx] * gainratios[gainidx] * pedrmss[gainidx] * pedrmss[gainidx] *
0155                       pedestal.asDiagonal() * noisecors[gainidx] * pedestal.asDiagonal();
0156           if (!dynamicPedestal && _addPedestalUncertainty > 0.) {
0157             //add fully correlated component to noise covariance to inflate pedestal uncertainty
0158             noisecov += gainratios[gainidx] * gainratios[gainidx] * _addPedestalUncertainty * _addPedestalUncertainty *
0159                         pedestal.asDiagonal() * SampleMatrix::Ones() * pedestal.asDiagonal();
0160           }
0161         }
0162       }
0163     }
0164   } else {
0165     noisecov = aped->rms_x12 * aped->rms_x12 * noisecors[0];
0166     if (!dynamicPedestal && _addPedestalUncertainty > 0.) {
0167       //add fully correlated component to noise covariance to inflate pedestal uncertainty
0168       noisecov += _addPedestalUncertainty * _addPedestalUncertainty * SampleMatrix::Ones();
0169     }
0170   }
0171
0172   //optimized one-pulse fit for hlt
0173   bool usePrefit = false;
0174   if (_doPrefit) {
0175     status =
0176         _pulsefuncSingle.DoFit(amplitudes, noisecov, _singlebx, fullpulse, fullpulsecov, gainsPedestal, badSamples);
0177     amplitude = status ? _pulsefuncSingle.X()[0] : 0.;
0178     amperr = status ? _pulsefuncSingle.Errors()[0] : 0.;
0179     chisq = _pulsefuncSingle.ChiSq();
0180
0181     if (chisq < _prefitMaxChiSq) {
0182       usePrefit = true;
0183     }
0184   }
0185
0186   if (!usePrefit) {
0187     if (!_computeErrors)
0188       _pulsefunc.disableErrorCalculation();
0189     status = _pulsefunc.DoFit(amplitudes, noisecov, activeBX, fullpulse, fullpulsecov, gainsPedestal, badSamples);
0190     chisq = _pulsefunc.ChiSq();
0191
0192     if (!status) {
0193       edm::LogWarning("EcalUncalibRecHitMultiFitAlgo::makeRecHit") << "Failed Fit" << std::endl;
0194     }
0195
0196     unsigned int ipulseintime = 0;
0197     for (unsigned int ipulse = 0; ipulse < _pulsefunc.BXs().rows(); ++ipulse) {
0198       if (_pulsefunc.BXs().coeff(ipulse) == 0) {
0199         ipulseintime = ipulse;
0200         break;
0201       }
0202     }
0203
0204     amplitude = status ? _pulsefunc.X()[ipulseintime] : 0.;
0205     amperr = status ? _pulsefunc.Errors()[ipulseintime] : 0.;
0206   }
0207
0208   double jitter = 0.;
0209
0210   EcalUncalibratedRecHit rh(dataFrame.id(), amplitude, pedval, jitter, chisq, flags);
0211   rh.setAmplitudeError(amperr);
0212
0213   if (!usePrefit) {
0214     for (unsigned int ipulse = 0; ipulse < _pulsefunc.BXs().rows(); ++ipulse) {
0215       int bx = _pulsefunc.BXs().coeff(ipulse);
0216       if (bx != 0 && std::abs(bx) < 100) {
0217         rh.setOutOfTimeAmplitude(bx + 5, status ? _pulsefunc.X().coeff(ipulse) : 0.);
0218       } else if (bx == (100 + gainsPedestal[iSampleMax])) {
0219         rh.setPedestal(status ? _pulsefunc.X().coeff(ipulse) : 0.);
0220       }
0221     }
0222   }
0223
0224   return rh;
0225 }