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