CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EcalUncalibRecHitWorkerGlobal.cc
Go to the documentation of this file.
1 #include "RecoLocalCalo/EcalRecProducers/plugins/EcalUncalibRecHitWorkerGlobal.h"
2 
3 #include "FWCore/Framework/interface/EventSetup.h"
4 #include "FWCore/Framework/interface/Event.h"
5 #include "FWCore/MessageLogger/interface/MessageLogger.h"
6 #include "FWCore/ParameterSet/interface/ParameterSet.h"
7 
8 #include "CondFormats/DataRecord/interface/EcalGainRatiosRcd.h"
9 #include "CondFormats/DataRecord/interface/EcalPedestalsRcd.h"
10 #include "CondFormats/DataRecord/interface/EcalWeightXtalGroupsRcd.h"
11 #include "CondFormats/DataRecord/interface/EcalTBWeightsRcd.h"
12 #include "CondFormats/DataRecord/interface/EcalSampleMaskRcd.h"
13 #include "CondFormats/DataRecord/interface/EcalTimeCalibConstantsRcd.h"
14 #include "CondFormats/DataRecord/interface/EcalTimeOffsetConstantRcd.h"
15 #include "CondFormats/DataRecord/interface/EcalTimeBiasCorrectionsRcd.h"
16 
17 #include <FWCore/ParameterSet/interface/ConfigurationDescriptions.h>
18 #include <FWCore/ParameterSet/interface/ParameterSetDescription.h>
19 
20 EcalUncalibRecHitWorkerGlobal::EcalUncalibRecHitWorkerGlobal(const edm::ParameterSet&ps,edm::ConsumesCollector& c) :
21  EcalUncalibRecHitWorkerBaseClass(ps,c)
22 {
23  // ratio method parameters
24  EBtimeFitParameters_ = ps.getParameter<std::vector<double> >("EBtimeFitParameters");
25  EEtimeFitParameters_ = ps.getParameter<std::vector<double> >("EEtimeFitParameters");
26  EBamplitudeFitParameters_ = ps.getParameter<std::vector<double> >("EBamplitudeFitParameters");
27  EEamplitudeFitParameters_ = ps.getParameter<std::vector<double> >("EEamplitudeFitParameters");
28  EBtimeFitLimits_.first = ps.getParameter<double>("EBtimeFitLimits_Lower");
29  EBtimeFitLimits_.second = ps.getParameter<double>("EBtimeFitLimits_Upper");
30  EEtimeFitLimits_.first = ps.getParameter<double>("EEtimeFitLimits_Lower");
31  EEtimeFitLimits_.second = ps.getParameter<double>("EEtimeFitLimits_Upper");
32  EBtimeConstantTerm_=ps.getParameter<double>("EBtimeConstantTerm");
33  EBtimeNconst_=ps.getParameter<double>("EBtimeNconst");
34  EEtimeConstantTerm_=ps.getParameter<double>("EEtimeConstantTerm");
35  EEtimeNconst_=ps.getParameter<double>("EEtimeNconst");
36  outOfTimeThreshG12pEB_ = ps.getParameter<double>("outOfTimeThresholdGain12pEB");
37  outOfTimeThreshG12mEB_ = ps.getParameter<double>("outOfTimeThresholdGain12mEB");
38  outOfTimeThreshG61pEB_ = ps.getParameter<double>("outOfTimeThresholdGain61pEB");
39  outOfTimeThreshG61mEB_ = ps.getParameter<double>("outOfTimeThresholdGain61mEB");
40  outOfTimeThreshG12pEE_ = ps.getParameter<double>("outOfTimeThresholdGain12pEE");
41  outOfTimeThreshG12mEE_ = ps.getParameter<double>("outOfTimeThresholdGain12mEE");
42  outOfTimeThreshG61pEE_ = ps.getParameter<double>("outOfTimeThresholdGain61pEE");
43  outOfTimeThreshG61mEE_ = ps.getParameter<double>("outOfTimeThresholdGain61mEE");
44  amplitudeThreshEB_ = ps.getParameter<double>("amplitudeThresholdEB");
45  amplitudeThreshEE_ = ps.getParameter<double>("amplitudeThresholdEE");
46 
47  // spike threshold
48  ebSpikeThresh_ = ps.getParameter<double>("ebSpikeThreshold");
49  // leading edge parameters
50  ebPulseShape_ = ps.getParameter<std::vector<double> >("ebPulseShape");
51  eePulseShape_ = ps.getParameter<std::vector<double> >("eePulseShape");
52  // chi2 parameters
53  kPoorRecoFlagEB_ = ps.getParameter<bool>("kPoorRecoFlagEB");
54  kPoorRecoFlagEE_ = ps.getParameter<bool>("kPoorRecoFlagEE");;
55  chi2ThreshEB_=ps.getParameter<double>("chi2ThreshEB_");
56  chi2ThreshEE_=ps.getParameter<double>("chi2ThreshEE_");
57  EBchi2Parameters_ = ps.getParameter<std::vector<double> >("EBchi2Parameters");
58  EEchi2Parameters_ = ps.getParameter<std::vector<double> >("EEchi2Parameters");
59 }
60 
61 
62 
63 EcalUncalibRecHitWorkerGlobal::EcalUncalibRecHitWorkerGlobal(const edm::ParameterSet&ps) :
64  EcalUncalibRecHitWorkerBaseClass(ps)
65 {
66  // ratio method parameters
67  EBtimeFitParameters_ = ps.getParameter<std::vector<double> >("EBtimeFitParameters");
68  EEtimeFitParameters_ = ps.getParameter<std::vector<double> >("EEtimeFitParameters");
69  EBamplitudeFitParameters_ = ps.getParameter<std::vector<double> >("EBamplitudeFitParameters");
70  EEamplitudeFitParameters_ = ps.getParameter<std::vector<double> >("EEamplitudeFitParameters");
71  EBtimeFitLimits_.first = ps.getParameter<double>("EBtimeFitLimits_Lower");
72  EBtimeFitLimits_.second = ps.getParameter<double>("EBtimeFitLimits_Upper");
73  EEtimeFitLimits_.first = ps.getParameter<double>("EEtimeFitLimits_Lower");
74  EEtimeFitLimits_.second = ps.getParameter<double>("EEtimeFitLimits_Upper");
75  EBtimeConstantTerm_=ps.getParameter<double>("EBtimeConstantTerm");
76  EBtimeNconst_=ps.getParameter<double>("EBtimeNconst");
77  EEtimeConstantTerm_=ps.getParameter<double>("EEtimeConstantTerm");
78  EEtimeNconst_=ps.getParameter<double>("EEtimeNconst");
79  outOfTimeThreshG12pEB_ = ps.getParameter<double>("outOfTimeThresholdGain12pEB");
80  outOfTimeThreshG12mEB_ = ps.getParameter<double>("outOfTimeThresholdGain12mEB");
81  outOfTimeThreshG61pEB_ = ps.getParameter<double>("outOfTimeThresholdGain61pEB");
82  outOfTimeThreshG61mEB_ = ps.getParameter<double>("outOfTimeThresholdGain61mEB");
83  outOfTimeThreshG12pEE_ = ps.getParameter<double>("outOfTimeThresholdGain12pEE");
84  outOfTimeThreshG12mEE_ = ps.getParameter<double>("outOfTimeThresholdGain12mEE");
85  outOfTimeThreshG61pEE_ = ps.getParameter<double>("outOfTimeThresholdGain61pEE");
86  outOfTimeThreshG61mEE_ = ps.getParameter<double>("outOfTimeThresholdGain61mEE");
87  amplitudeThreshEB_ = ps.getParameter<double>("amplitudeThresholdEB");
88  amplitudeThreshEE_ = ps.getParameter<double>("amplitudeThresholdEE");
89 
90  // spike threshold
91  ebSpikeThresh_ = ps.getParameter<double>("ebSpikeThreshold");
92  // leading edge parameters
93  ebPulseShape_ = ps.getParameter<std::vector<double> >("ebPulseShape");
94  eePulseShape_ = ps.getParameter<std::vector<double> >("eePulseShape");
95  // chi2 parameters
96  kPoorRecoFlagEB_ = ps.getParameter<bool>("kPoorRecoFlagEB");
97  kPoorRecoFlagEE_ = ps.getParameter<bool>("kPoorRecoFlagEE");;
98  chi2ThreshEB_=ps.getParameter<double>("chi2ThreshEB_");
99  chi2ThreshEE_=ps.getParameter<double>("chi2ThreshEE_");
100  EBchi2Parameters_ = ps.getParameter<std::vector<double> >("EBchi2Parameters");
101  EEchi2Parameters_ = ps.getParameter<std::vector<double> >("EEchi2Parameters");
102 }
103 
104 void
105 EcalUncalibRecHitWorkerGlobal::set(const edm::EventSetup& es)
106 {
107  // common setup
108  es.get<EcalGainRatiosRcd>().get(gains);
109  es.get<EcalPedestalsRcd>().get(peds);
110 
111  // for the weights method
112  es.get<EcalWeightXtalGroupsRcd>().get(grps);
113  es.get<EcalTBWeightsRcd>().get(wgts);
114 
115  // which of the samples need be used
116  es.get<EcalSampleMaskRcd>().get(sampleMaskHand_);
117 
118  // for the ratio method
119 
120  // for the leading edge method
121  es.get<EcalTimeCalibConstantsRcd>().get(itime);
122  es.get<EcalTimeOffsetConstantRcd>().get(offtime);
123 
124  // for the time correction methods
125  es.get<EcalTimeBiasCorrectionsRcd>().get(timeCorrBias_);
126 }
127 
128 
129 // check saturation: 5 samples with gainId = 0
130 template < class C >
131 int EcalUncalibRecHitWorkerGlobal::isSaturated(const C & dataFrame)
132 {
133  //bool saturated_ = 0;
134  int cnt;
135  for (int j = 0; j < C::MAXSAMPLES - 5; ++j) {
136  cnt = 0;
137  for (int i = j; i < (j + 5) && i < C::MAXSAMPLES; ++i) {
138  if ( dataFrame.sample(i).gainId() == 0 ) ++cnt;
139  }
140  if ( cnt == 5 ) return j-1 ; // the last unsaturated sample
141  }
142  return -1; // no saturation found
143 }
144 
154 double EcalUncalibRecHitWorkerGlobal::timeCorrection(
155  float ampli,
156  const std::vector<float>& amplitudeBins,
157  const std::vector<float>& shiftBins) {
158 
159  // computed initially in ns. Than turned in the BX's, as
160  // EcalUncalibratedRecHit need be.
161  double theCorrection = 0;
162 
163  // sanity check for arrays
164  if (amplitudeBins.size() == 0) {
165  edm::LogError("EcalRecHitError")
166  << "timeCorrAmplitudeBins is empty, forcing no time bias corrections.";
167 
168  return 0;
169  }
170 
171  if (amplitudeBins.size() != shiftBins.size()) {
172  edm::LogError("EcalRecHitError")
173  << "Size of timeCorrAmplitudeBins different from "
174  "timeCorrShiftBins. Forcing no time bias corrections. ";
175 
176  return 0;
177  }
178 
179  int myBin = -1;
180  for (int bin = 0; bin < (int) amplitudeBins.size(); bin++) {
181  if (ampli > amplitudeBins.at(bin)) {
182  myBin = bin;
183  } else {
184  break;
185  }
186  }
187 
188  if (myBin == -1) {
189  theCorrection = shiftBins.at(0);
190  } else if (myBin == ((int)(amplitudeBins.size() - 1))) {
191  theCorrection = shiftBins.at(myBin);
192  } else if (-1 < myBin && myBin < ((int) amplitudeBins.size() - 1)) {
193  // interpolate linearly between two assingned points
194  theCorrection = (shiftBins.at(myBin + 1) - shiftBins.at(myBin));
195  theCorrection *= (((double) ampli) - amplitudeBins.at(myBin)) /
196  (amplitudeBins.at(myBin + 1) - amplitudeBins.at(myBin));
197  theCorrection += shiftBins.at(myBin);
198  } else {
199  edm::LogError("EcalRecHitError")
200  << "Assigning time correction impossible. Setting it to 0 ";
201  theCorrection = 0.;
202  }
203 
204  // convert ns into clocks
205  return theCorrection / 25.;
206 }
207 
208 
209 
210 bool
211 EcalUncalibRecHitWorkerGlobal::run( const edm::Event & evt,
212  const EcalDigiCollection::const_iterator & itdg,
213  EcalUncalibratedRecHitCollection & result )
214 {
215  DetId detid(itdg->id());
216 
217  const EcalSampleMask *sampleMask_ = sampleMaskHand_.product();
218 
219  // intelligence for recHit computation
220  EcalUncalibratedRecHit uncalibRecHit;
221 
222 
223  const EcalPedestals::Item * aped = 0;
224  const EcalMGPAGainRatio * aGain = 0;
225  const EcalXtalGroupId * gid = 0;
226  float offsetTime = 0;
227 
228  if (detid.subdetId()==EcalEndcap) {
229  unsigned int hashedIndex = EEDetId(detid).hashedIndex();
230  aped = &peds->endcap(hashedIndex);
231  aGain = &gains->endcap(hashedIndex);
232  gid = &grps->endcap(hashedIndex);
233  offsetTime = offtime->getEEValue();
234  } else {
235  unsigned int hashedIndex = EBDetId(detid).hashedIndex();
236  aped = &peds->barrel(hashedIndex);
237  aGain = &gains->barrel(hashedIndex);
238  gid = &grps->barrel(hashedIndex);
239  offsetTime = offtime->getEBValue();
240  }
241 
242  pedVec[0] = aped->mean_x12;
243  pedVec[1] = aped->mean_x6;
244  pedVec[2] = aped->mean_x1;
245  pedRMSVec[0] = aped->rms_x12;
246  pedRMSVec[1] = aped->rms_x6;
247  pedRMSVec[2] = aped->rms_x1;
248  gainRatios[0] = 1.;
249  gainRatios[1] = aGain->gain12Over6();
250  gainRatios[2] = aGain->gain6Over1()*aGain->gain12Over6();
251 
252  // compute the right bin of the pulse shape using time calibration constants
253  EcalTimeCalibConstantMap::const_iterator it = itime->find( detid );
254  EcalTimeCalibConstant itimeconst = 0;
255  if( it != itime->end() ) {
256  itimeconst = (*it);
257  } else {
258  edm::LogError("EcalRecHitError") << "No time intercalib const found for xtal "
259  << detid.rawId()
260  << "! something wrong with EcalTimeCalibConstants in your DB? ";
261  }
262 
263 
264  // === amplitude computation ===
265  int leadingSample = -1;
266  if (detid.subdetId()==EcalEndcap) {
267  leadingSample = ((EcalDataFrame)(*itdg)).lastUnsaturatedSample();
268  } else {
269  leadingSample = ((EcalDataFrame)(*itdg)).lastUnsaturatedSample();
270  }
271 
272  if ( leadingSample >= 0 ) { // saturation
273  if ( leadingSample != 4 ) {
274  // all samples different from the fifth are not reliable for the amplitude estimation
275  // put by default the energy at the saturation threshold and flag as saturated
276  float sratio = 1;
277  if ( detid.subdetId()==EcalBarrel) {
278  sratio = ebPulseShape_[5] / ebPulseShape_[4];
279  } else {
280  sratio = eePulseShape_[5] / eePulseShape_[4];
281  }
282  uncalibRecHit = EcalUncalibratedRecHit( (*itdg).id(), 4095*12*sratio, 0, 0, 0);
283  uncalibRecHit.setFlagBit( EcalUncalibratedRecHit::kSaturated );
284  } else {
285  // float clockToNsConstant = 25.;
286  // reconstruct the rechit
287  if (detid.subdetId()==EcalEndcap) {
288  leadingEdgeMethod_endcap_.setPulseShape( eePulseShape_ );
289  // float mult = (float)eePulseShape_.size() / (float)(*itdg).size();
290  // bin (or some analogous mapping) will be used instead of the leadingSample
291  //int bin = (int)(( (mult * leadingSample + mult/2) * clockToNsConstant + itimeconst ) / clockToNsConstant);
292  // bin is not uset for the moment
293  leadingEdgeMethod_endcap_.setLeadingEdgeSample( leadingSample );
294  uncalibRecHit = leadingEdgeMethod_endcap_.makeRecHit(*itdg, pedVec, gainRatios, 0, 0);
295  uncalibRecHit.setFlagBit( EcalUncalibratedRecHit::kLeadingEdgeRecovered );
296  leadingEdgeMethod_endcap_.setLeadingEdgeSample( -1 );
297  } else {
298  leadingEdgeMethod_barrel_.setPulseShape( ebPulseShape_ );
299  // float mult = (float)ebPulseShape_.size() / (float)(*itdg).size();
300  // bin (or some analogous mapping) will be used instead of the leadingSample
301  //int bin = (int)(( (mult * leadingSample + mult/2) * clockToNsConstant + itimeconst ) / clockToNsConstant);
302  // bin is not uset for the moment
303  leadingEdgeMethod_barrel_.setLeadingEdgeSample( leadingSample );
304  uncalibRecHit = leadingEdgeMethod_barrel_.makeRecHit(*itdg, pedVec, gainRatios, 0, 0);
305  uncalibRecHit.setFlagBit( EcalUncalibratedRecHit::kLeadingEdgeRecovered );
306  leadingEdgeMethod_barrel_.setLeadingEdgeSample( -1 );
307  }
308  }
309  // do not propagate the default chi2 = -1 value to the calib rechit (mapped to 64), set it to 0 when saturation
310  uncalibRecHit.setChi2(0);
311  } else {
312  // weights method
313  EcalTBWeights::EcalTDCId tdcid(1);
314  EcalTBWeights::EcalTBWeightMap const & wgtsMap = wgts->getMap();
315  EcalTBWeights::EcalTBWeightMap::const_iterator wit;
316  wit = wgtsMap.find( std::make_pair(*gid,tdcid) );
317  if( wit == wgtsMap.end() ) {
318  edm::LogError("EcalUncalibRecHitError") << "No weights found for EcalGroupId: "
319  << gid->id() << " and EcalTDCId: " << tdcid
320  << "\n skipping digi with id: " << detid.rawId();
321 
322  return false;
323  }
324  const EcalWeightSet& wset = wit->second; // this is the EcalWeightSet
325 
326  const EcalWeightSet::EcalWeightMatrix& mat1 = wset.getWeightsBeforeGainSwitch();
327  const EcalWeightSet::EcalWeightMatrix& mat2 = wset.getWeightsAfterGainSwitch();
328 
329  weights[0] = &mat1;
330  weights[1] = &mat2;
331 
332  // get uncalibrated recHit from weights
333  if (detid.subdetId()==EcalEndcap) {
334  uncalibRecHit = weightsMethod_endcap_.makeRecHit(*itdg, pedVec, pedRMSVec, gainRatios, weights, testbeamEEShape);
335  } else {
336  uncalibRecHit = weightsMethod_barrel_.makeRecHit(*itdg, pedVec, pedRMSVec, gainRatios, weights, testbeamEBShape);
337  }
338 
339  // === time computation ===
340  // ratio method
341  float const clockToNsConstant = 25.;
342  if (detid.subdetId()==EcalEndcap) {
343  ratioMethod_endcap_.init( *itdg, *sampleMask_, pedVec, pedRMSVec, gainRatios );
344  ratioMethod_endcap_.computeTime( EEtimeFitParameters_, EEtimeFitLimits_, EEamplitudeFitParameters_ );
345  ratioMethod_endcap_.computeAmplitude( EEamplitudeFitParameters_);
346  EcalUncalibRecHitRatioMethodAlgo<EEDataFrame>::CalculatedRecHit crh = ratioMethod_endcap_.getCalculatedRecHit();
347  double theTimeCorrectionEE = timeCorrection(uncalibRecHit.amplitude(),
348  timeCorrBias_->EETimeCorrAmplitudeBins, timeCorrBias_->EETimeCorrShiftBins);
349 
350  uncalibRecHit.setJitter( crh.timeMax - 5 + theTimeCorrectionEE);
351  uncalibRecHit.setJitterError( std::sqrt(pow(crh.timeError,2) + std::pow(EEtimeConstantTerm_,2)/std::pow(clockToNsConstant,2)) );
352  // consider flagging as kOutOfTime only if above noise
353  if (uncalibRecHit.amplitude() > pedRMSVec[0] * amplitudeThreshEE_){
354  float outOfTimeThreshP = outOfTimeThreshG12pEE_;
355  float outOfTimeThreshM = outOfTimeThreshG12mEE_;
356  // determine if gain has switched away from gainId==1 (x12 gain)
357  // and determine cuts (number of 'sigmas') to ose for kOutOfTime
358  // >3k ADC is necessasry condition for gain switch to occur
359  if (uncalibRecHit.amplitude() > 3000.){
360  for (int iSample = 0; iSample < EEDataFrame::MAXSAMPLES; iSample++) {
361  int GainId = ((EcalDataFrame)(*itdg)).sample(iSample).gainId();
362  if (GainId!=1) {
363  outOfTimeThreshP = outOfTimeThreshG61pEE_;
364  outOfTimeThreshM = outOfTimeThreshG61mEE_;
365  break;
366  }
367  }}
368  float correctedTime = (crh.timeMax-5) * clockToNsConstant + itimeconst + offsetTime;
369  float cterm = EEtimeConstantTerm_;
370  float sigmaped = pedRMSVec[0]; // approx for lower gains
371  float nterm = EEtimeNconst_*sigmaped/uncalibRecHit.amplitude();
372  float sigmat = std::sqrt( nterm*nterm + cterm*cterm );
373  if ( ( correctedTime > sigmat*outOfTimeThreshP ) ||
374  ( correctedTime < (-1.*sigmat*outOfTimeThreshM) ))
375  { uncalibRecHit.setFlagBit( EcalUncalibratedRecHit::kOutOfTime ); }
376  }
377 
378  } else {
379  ratioMethod_barrel_.init( *itdg, *sampleMask_, pedVec, pedRMSVec, gainRatios );
380  ratioMethod_barrel_.fixMGPAslew(*itdg);
381  ratioMethod_barrel_.computeTime( EBtimeFitParameters_, EBtimeFitLimits_, EBamplitudeFitParameters_ );
382  ratioMethod_barrel_.computeAmplitude( EBamplitudeFitParameters_);
383  EcalUncalibRecHitRatioMethodAlgo<EBDataFrame>::CalculatedRecHit crh = ratioMethod_barrel_.getCalculatedRecHit();
384 
385  double theTimeCorrectionEB = timeCorrection(uncalibRecHit.amplitude(),
386  timeCorrBias_->EBTimeCorrAmplitudeBins, timeCorrBias_->EBTimeCorrShiftBins);
387 
388  uncalibRecHit.setJitter( crh.timeMax - 5 + theTimeCorrectionEB);
389 
390  uncalibRecHit.setJitterError( std::sqrt(std::pow(crh.timeError,2) + std::pow(EBtimeConstantTerm_,2)/std::pow(clockToNsConstant,2)) );
391  // consider flagging as kOutOfTime only if above noise
392  if (uncalibRecHit.amplitude() > pedRMSVec[0] * amplitudeThreshEB_){
393  float outOfTimeThreshP = outOfTimeThreshG12pEB_;
394  float outOfTimeThreshM = outOfTimeThreshG12mEB_;
395  // determine if gain has switched away from gainId==1 (x12 gain)
396  // and determine cuts (number of 'sigmas') to ose for kOutOfTime
397  // >3k ADC is necessasry condition for gain switch to occur
398  if (uncalibRecHit.amplitude() > 3000.){
399  for (int iSample = 0; iSample < EBDataFrame::MAXSAMPLES; iSample++) {
400  int GainId = ((EcalDataFrame)(*itdg)).sample(iSample).gainId();
401  if (GainId!=1) {
402  outOfTimeThreshP = outOfTimeThreshG61pEB_;
403  outOfTimeThreshM = outOfTimeThreshG61mEB_;
404  break;}
405  } }
406  float correctedTime = (crh.timeMax-5) * clockToNsConstant + itimeconst + offsetTime;
407  float cterm = EBtimeConstantTerm_;
408  float sigmaped = pedRMSVec[0]; // approx for lower gains
409  float nterm = EBtimeNconst_*sigmaped/uncalibRecHit.amplitude();
410  float sigmat = std::sqrt( nterm*nterm + cterm*cterm );
411  if ( ( correctedTime > sigmat*outOfTimeThreshP ) ||
412  ( correctedTime < (-1.*sigmat*outOfTimeThreshM) ))
413  { uncalibRecHit.setFlagBit( EcalUncalibratedRecHit::kOutOfTime ); }
414  }
415  }
416 
417  // === chi2express ===
418  if (detid.subdetId()==EcalEndcap) {
419 
420  double amplitude = uncalibRecHit.amplitude();
421  double amplitudeOutOfTime = 0.;
422  double jitter= uncalibRecHit.jitter();
423 
424 
425 
426  EcalUncalibRecHitRecChi2Algo<EEDataFrame>chi2expressEE_(
427  *itdg,
428  amplitude,
429  (itimeconst + offsetTime),
430  amplitudeOutOfTime,
431  jitter,
432  pedVec,
433  pedRMSVec,
434  gainRatios,
435  testbeamEEShape,
436  EEchi2Parameters_
437  );
438  double chi2 = chi2expressEE_.chi2();
439  uncalibRecHit.setChi2(chi2);
440 
441  if(kPoorRecoFlagEE_)
442  {
443 
444  if (chi2>chi2ThreshEE_) {
445 
446  // first check if all samples are ok, if not don't use chi2 to flag
447  bool samplesok = true;
448  for (int sample =0; sample < EcalDataFrame::MAXSAMPLES; ++sample) {
449  if (!sampleMask_->useSampleEE(sample)) {
450  samplesok = false;
451  break;
452  }
453  }
454  if (samplesok) uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kPoorReco);
455  }
456 
457 
458  }
459 
460  } else {
461  double amplitude = uncalibRecHit.amplitude();
462  double amplitudeOutOfTime = 0.;
463  double jitter= uncalibRecHit.jitter();
464 
465  EcalUncalibRecHitRecChi2Algo<EBDataFrame>chi2expressEB_(
466  *itdg,
467  amplitude,
468  (itimeconst + offsetTime),
469  amplitudeOutOfTime,
470  jitter,
471  pedVec,
472  pedRMSVec,
473  gainRatios,
474  testbeamEBShape,
475  EBchi2Parameters_
476  );
477  double chi2 = chi2expressEB_.chi2();
478  uncalibRecHit.setChi2(chi2);
479 
480  if(kPoorRecoFlagEB_)
481  {
482 
483  if(chi2>chi2ThreshEB_){
484  // first check if all samples are ok, if not don't use chi2 to flag
485  bool samplesok = true;
486  for (int sample =0; sample < EcalDataFrame::MAXSAMPLES; ++sample) {
487  if (!sampleMask_->useSampleEB(sample)) {
488  samplesok = false;
489  break;
490  }
491  }
492  if (samplesok) uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kPoorReco);
493  }
494 
495  }
496  }
497  }
498 
499  // set flags if gain switch has occurred
500  if( ((EcalDataFrame)(*itdg)).hasSwitchToGain6() ) uncalibRecHit.setFlagBit( EcalUncalibratedRecHit::kHasSwitchToGain6 );
501  if( ((EcalDataFrame)(*itdg)).hasSwitchToGain1() ) uncalibRecHit.setFlagBit( EcalUncalibratedRecHit::kHasSwitchToGain1 );
502 
503  // put the recHit in the collection
504  if (detid.subdetId()==EcalEndcap) {
505  result.push_back( uncalibRecHit );
506  } else {
507  result.push_back( uncalibRecHit );
508  }
509 
510  return true;
511 }
512 
513 edm::ParameterSetDescription
514 EcalUncalibRecHitWorkerGlobal::getAlgoDescription() {
515 
516  edm::ParameterSetDescription psd;
517  psd.addNode(edm::ParameterDescription<std::vector<double>>("eePulseShape", {5.2e-05,-5.26e-05,6.66e-05,0.1168,0.7575,1.0,0.8876,0.6732,0.4741,0.3194}, true) and
518  edm::ParameterDescription<std::vector<double>>("EBtimeFitParameters", {-2.015452,3.130702,-12.3473,41.88921,-82.83944,91.01147,-50.35761,11.05621}, true) and
519  edm::ParameterDescription<double>("outOfTimeThresholdGain61pEB", 5, true) and
520  edm::ParameterDescription<double>("amplitudeThresholdEE", 10, true) and
521  edm::ParameterDescription<double>("EBtimeConstantTerm", 0.6, true) and
522  edm::ParameterDescription<double>("outOfTimeThresholdGain61pEE", 1000, true) and
523  edm::ParameterDescription<double>("ebSpikeThreshold", 1.042, true) and
524  edm::ParameterDescription<double>("EBtimeNconst", 28.5, true) and
525  edm::ParameterDescription<bool>("kPoorRecoFlagEB", true, true) and
526  edm::ParameterDescription<std::vector<double>>("ebPulseShape", {5.2e-05,-5.26e-05,6.66e-05,0.1168,0.7575,1.0,0.8876,0.6732,0.4741,0.3194}, true) and
527  edm::ParameterDescription<double>("EBtimeFitLimits_Lower", 0.2, true) and
528  edm::ParameterDescription<bool>("kPoorRecoFlagEE", false, true) and
529  edm::ParameterDescription<double>("chi2ThreshEB_", 36.0, true) and
530  edm::ParameterDescription<std::vector<double>>("EEtimeFitParameters", {-2.390548,3.553628,-17.62341,67.67538,-133.213,140.7432,-75.41106,16.20277}, true) and
531  edm::ParameterDescription<double>("outOfTimeThresholdGain61mEE", 1000, true) and
532  edm::ParameterDescription<std::vector<double>>("EEchi2Parameters", {2.122,0.022,2.122,0.022}, true) and
533  edm::ParameterDescription<double>("outOfTimeThresholdGain12mEE", 1000, true) and
534  edm::ParameterDescription<double>("outOfTimeThresholdGain12mEB", 5, true) and
535  edm::ParameterDescription<double>("EEtimeFitLimits_Upper", 1.4, true) and
536  edm::ParameterDescription<double>("EEtimeFitLimits_Lower", 0.2, true) and
537  edm::ParameterDescription<std::vector<double>>("EEamplitudeFitParameters", {1.89,1.4}, true) and
538  edm::ParameterDescription<std::vector<double>>("EBamplitudeFitParameters", {1.138,1.652}, true) and
539  edm::ParameterDescription<double>("amplitudeThresholdEB", 10, true) and
540  edm::ParameterDescription<double>("outOfTimeThresholdGain12pEE", 1000, true) and
541  edm::ParameterDescription<double>("outOfTimeThresholdGain12pEB", 5, true) and
542  edm::ParameterDescription<double>("EEtimeNconst", 31.8, true) and
543  edm::ParameterDescription<double>("outOfTimeThresholdGain61mEB", 5, true) and
544  edm::ParameterDescription<std::vector<double>>("EBchi2Parameters", {2.122,0.022,2.122,0.022}, true) and
545  edm::ParameterDescription<double>("EEtimeConstantTerm", 1.0, true) and
546  edm::ParameterDescription<double>("chi2ThreshEE_", 95.0, true) and
547  edm::ParameterDescription<double>("EBtimeFitLimits_Upper", 1.4, true));
548 
549  return psd;
550 }
551 
552 
553 #include "FWCore/Framework/interface/MakerMacros.h"
554 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitWorkerFactory.h"
555 DEFINE_EDM_PLUGIN( EcalUncalibRecHitWorkerFactory, EcalUncalibRecHitWorkerGlobal, "EcalUncalibRecHitWorkerGlobal" );
556 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitFillDescriptionWorkerFactory.h"
557 DEFINE_EDM_PLUGIN( EcalUncalibRecHitFillDescriptionWorkerFactory, EcalUncalibRecHitWorkerGlobal, "EcalUncalibRecHitWorkerGlobal");
EcalUncalibRecHitWorkerGlobal::offtime
edm::ESHandle< EcalTimeOffsetConstant > offtime
Definition: EcalUncalibRecHitWorkerGlobal.h:107
EcalTBWeights::EcalTBWeightMap
std::map< std::pair< EcalXtalGroupId, EcalTDCId >, EcalWeightSet > EcalTBWeightMap
Definition: EcalTBWeights.h:21
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
i
int i
Definition: DBlmapReader.cc:9
EcalUncalibRecHitWorkerGlobal::timeCorrBias_
edm::ESHandle< EcalTimeBiasCorrections > timeCorrBias_
Definition: EcalUncalibRecHitWorkerGlobal.h:103
EcalUncalibRecHitWorkerGlobal::run
bool run(const edm::Event &evt, const EcalDigiCollection::const_iterator &digi, EcalUncalibratedRecHitCollection &result)
Definition: EcalUncalibRecHitWorkerGlobal.cc:211
EBDetId::hashedIndex
int hashedIndex() const
get a compact index for arrays
Definition: EBDetId.h:86
EcalUncalibRecHitRatioMethodAlgo::computeAmplitude
void computeAmplitude(std::vector< double > &amplitudeFitParameters)
Definition: EcalUncalibRecHitRatioMethodAlgo.h:560
edm::DataFrameContainer::const_iterator
boost::transform_iterator< IterHelp, boost::counting_iterator< int > > const_iterator
Definition: DataFrameContainer.h:51
EcalUncalibRecHitWorkerGlobal::timeCorrection
double timeCorrection(float ampli, const std::vector< float > &amplitudeBins, const std::vector< float > &shiftBins)
Definition: EcalUncalibRecHitWorkerGlobal.cc:154
EcalUncalibratedRecHit::setJitter
void setJitter(float jitter)
Definition: EcalUncalibratedRecHit.h:46
EcalUncalibRecHitWorkerGlobal::EBtimeFitLimits_
std::pair< double, double > EBtimeFitLimits_
Definition: EcalUncalibRecHitWorkerGlobal.h:81
EcalUncalibRecHitRatioMethodAlgo::computeTime
void computeTime(std::vector< double > &timeFitParameters, std::pair< double, double > &timeFitLimits, std::vector< double > &amplitudeFitParameters)
Definition: EcalUncalibRecHitRatioMethodAlgo.h:224
edm::SortedCollection::push_back
void push_back(T const &t)
Definition: SortedCollection.h:199
EcalUncalibRecHitWorkerGlobal::weightsMethod_barrel_
EcalUncalibRecHitRecWeightsAlgo< EBDataFrame > weightsMethod_barrel_
Definition: EcalUncalibRecHitWorkerGlobal.h:68
edm::ParameterSetDescription::addNode
ParameterDescriptionNode * addNode(ParameterDescriptionNode const &node)
Definition: ParameterSetDescription.cc:59
EcalUncalibratedRecHit::setJitterError
void setJitterError(float jitterErr)
Definition: EcalUncalibratedRecHit.cc:41
EcalUncalibRecHitWorkerGlobal::EEtimeFitLimits_
std::pair< double, double > EEtimeFitLimits_
Definition: EcalUncalibRecHitWorkerGlobal.h:82
EcalUncalibRecHitLeadingEdgeAlgo::setPulseShape
void setPulseShape(std::vector< double > &shape)
Definition: EcalUncalibRecHitLeadingEdgeAlgo.h:31
EcalUncalibRecHitWorkerGlobal::ratioMethod_endcap_
EcalUncalibRecHitRatioMethodAlgo< EEDataFrame > ratioMethod_endcap_
Definition: EcalUncalibRecHitWorkerGlobal.h:85
EcalXtalGroupId::id
const unsigned int id() const
Definition: EcalXtalGroupId.h:22
EcalUncalibRecHitWorkerGlobal::leadingEdgeMethod_barrel_
EcalUncalibRecHitLeadingEdgeAlgo< EBDataFrame > leadingEdgeMethod_barrel_
Definition: EcalUncalibRecHitWorkerGlobal.h:110
ecalpyutils::hashedIndex
int hashedIndex(int ieta, int iphi)
Definition: EcalPyUtils.cc:42
EcalUncalibRecHitWorkerGlobal::ratioMethod_barrel_
EcalUncalibRecHitRatioMethodAlgo< EBDataFrame > ratioMethod_barrel_
Definition: EcalUncalibRecHitWorkerGlobal.h:84
EcalUncalibRecHitWorkerGlobal::gains
edm::ESHandle< EcalGainRatios > gains
Definition: EcalUncalibRecHitWorkerGlobal.h:56
EcalUncalibRecHitWorkerGlobal::leadingEdgeMethod_endcap_
EcalUncalibRecHitLeadingEdgeAlgo< EEDataFrame > leadingEdgeMethod_endcap_
Definition: EcalUncalibRecHitWorkerGlobal.h:111
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
EcalUncalibRecHitWorkerGlobal::peds
edm::ESHandle< EcalPedestals > peds
Definition: EcalUncalibRecHitWorkerGlobal.h:55
query.result
tuple result
Definition: query.py:137
EcalUncalibRecHitRecWeightsAlgo::makeRecHit
virtual EcalUncalibratedRecHit makeRecHit(const C &dataFrame, const double *pedestals, const double *pedestalsRMS, const double *gainRatios, const EcalWeightSet::EcalWeightMatrix **weights, const EcalShapeBase &testbeamPulseShape)
Compute parameters.
Definition: EcalUncalibRecHitRecWeightsAlgo.h:30
j
int j
Definition: DBlmapReader.cc:9
EcalUncalibRecHitWorkerGlobal::wgts
edm::ESHandle< EcalTBWeights > wgts
Definition: EcalUncalibRecHitWorkerGlobal.h:65
EcalMGPAGainRatio::gain6Over1
float gain6Over1() const
Definition: EcalMGPAGainRatio.h:21
EcalUncalibRecHitWorkerGlobal::getAlgoDescription
edm::ParameterSetDescription getAlgoDescription()
Definition: EcalUncalibRecHitWorkerGlobal.cc:514
EcalUncalibRecHitWorkerGlobal::grps
edm::ESHandle< EcalWeightXtalGroups > grps
Definition: EcalUncalibRecHitWorkerGlobal.h:64
EcalWeightSet::getWeightsAfterGainSwitch
EcalWeightMatrix & getWeightsAfterGainSwitch()
Definition: EcalWeightSet.h:30
EcalUncalibRecHitLeadingEdgeAlgo::makeRecHit
virtual EcalUncalibratedRecHit makeRecHit(const C &dataFrame, const double *pedestals, const double *gainRatios, const EcalWeightSet::EcalWeightMatrix **weights, const EcalWeightSet::EcalChi2WeightMatrix **chi2Matrix)
Compute parameters.
Definition: EcalUncalibRecHitLeadingEdgeAlgo.h:35
DetId
Definition: DetId.h:18
EEDetId::hashedIndex
int hashedIndex() const
Definition: EEDetId.h:182
EcalWeightSet::getWeightsBeforeGainSwitch
EcalWeightMatrix & getWeightsBeforeGainSwitch()
Definition: EcalWeightSet.h:29
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
EcalCondObjectContainer< float >::const_iterator
std::vector< Item >::const_iterator const_iterator
Definition: EcalCondObjectContainer.h:17
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
EcalMGPAGainRatio::gain12Over6
float gain12Over6() const
Definition: EcalMGPAGainRatio.h:20
EcalUncalibRecHitWorkerGlobal::weights
const EcalWeightSet::EcalWeightMatrix * weights[2]
Definition: EcalUncalibRecHitWorkerGlobal.h:66
EcalTimeCalibConstant
float EcalTimeCalibConstant
Definition: EcalTimeCalibConstants.h:10
EcalUncalibRecHitWorkerGlobal::set
void set(const edm::EventSetup &es)
Definition: EcalUncalibRecHitWorkerGlobal.cc:105
EcalUncalibRecHitWorkerGlobal::itime
edm::ESHandle< EcalTimeCalibConstants > itime
Definition: EcalUncalibRecHitWorkerGlobal.h:106
EcalUncalibRecHitWorkerGlobal::weightsMethod_endcap_
EcalUncalibRecHitRecWeightsAlgo< EEDataFrame > weightsMethod_endcap_
Definition: EcalUncalibRecHitWorkerGlobal.h:69
EcalWeightSet::EcalWeightMatrix
math::Matrix< 3, 10 >::type EcalWeightMatrix
Definition: EcalWeightSet.h:22
DEFINE_EDM_PLUGIN
#define DEFINE_EDM_PLUGIN(factory, type, name)
Definition: PluginFactory.h:101
EcalUncalibRecHitWorkerGlobal::sampleMaskHand_
edm::ESHandle< EcalSampleMask > sampleMaskHand_
Definition: EcalUncalibRecHitWorkerGlobal.h:74
EcalUncalibRecHitRatioMethodAlgo::init
void init(const C &dataFrame, const EcalSampleMask &sampleMask, const double *pedestals, const double *pedestalRMSes, const double *gainRatios)
Definition: EcalUncalibRecHitRatioMethodAlgo.h:95
EcalDataFrame::MAXSAMPLES
static const int MAXSAMPLES
Definition: EcalDataFrame.h:48
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
compare_using_db.sample
string sample
Definition: compare_using_db.py:26
Generated for CMSSW Reference Manual by doxygen 1.8.5