00001 #include "RecoLocalCalo/EcalRecProducers/plugins/EcalUncalibRecHitWorkerGlobal.h"
00002
00003 #include "FWCore/Framework/interface/EventSetup.h"
00004 #include "FWCore/Framework/interface/Event.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00007
00008 #include "CondFormats/DataRecord/interface/EcalGainRatiosRcd.h"
00009 #include "CondFormats/DataRecord/interface/EcalPedestalsRcd.h"
00010 #include "CondFormats/DataRecord/interface/EcalWeightXtalGroupsRcd.h"
00011 #include "CondFormats/DataRecord/interface/EcalTBWeightsRcd.h"
00012 #include "CondFormats/DataRecord/interface/EcalTimeCalibConstantsRcd.h"
00013 #include "CondFormats/DataRecord/interface/EcalTimeOffsetConstantRcd.h"
00014
00015 EcalUncalibRecHitWorkerGlobal::EcalUncalibRecHitWorkerGlobal(const edm::ParameterSet&ps) :
00016 EcalUncalibRecHitWorkerBaseClass(ps)
00017 {
00018
00019 EBtimeFitParameters_ = ps.getParameter<std::vector<double> >("EBtimeFitParameters");
00020 EEtimeFitParameters_ = ps.getParameter<std::vector<double> >("EEtimeFitParameters");
00021 EBamplitudeFitParameters_ = ps.getParameter<std::vector<double> >("EBamplitudeFitParameters");
00022 EEamplitudeFitParameters_ = ps.getParameter<std::vector<double> >("EEamplitudeFitParameters");
00023 EBtimeFitLimits_.first = ps.getParameter<double>("EBtimeFitLimits_Lower");
00024 EBtimeFitLimits_.second = ps.getParameter<double>("EBtimeFitLimits_Upper");
00025 EEtimeFitLimits_.first = ps.getParameter<double>("EEtimeFitLimits_Lower");
00026 EEtimeFitLimits_.second = ps.getParameter<double>("EEtimeFitLimits_Upper");
00027 EBtimeConstantTerm_=ps.getParameter<double>("EBtimeConstantTerm");
00028 EBtimeNconst_=ps.getParameter<double>("EBtimeNconst");
00029 EEtimeConstantTerm_=ps.getParameter<double>("EEtimeConstantTerm");
00030 EEtimeNconst_=ps.getParameter<double>("EEtimeNconst");
00031 outOfTimeThreshG12pEB_ = ps.getParameter<double>("outOfTimeThresholdGain12pEB");
00032 outOfTimeThreshG12mEB_ = ps.getParameter<double>("outOfTimeThresholdGain12mEB");
00033 outOfTimeThreshG61pEB_ = ps.getParameter<double>("outOfTimeThresholdGain61pEB");
00034 outOfTimeThreshG61mEB_ = ps.getParameter<double>("outOfTimeThresholdGain61mEB");
00035 outOfTimeThreshG12pEE_ = ps.getParameter<double>("outOfTimeThresholdGain12pEE");
00036 outOfTimeThreshG12mEE_ = ps.getParameter<double>("outOfTimeThresholdGain12mEE");
00037 outOfTimeThreshG61pEE_ = ps.getParameter<double>("outOfTimeThresholdGain61pEE");
00038 outOfTimeThreshG61mEE_ = ps.getParameter<double>("outOfTimeThresholdGain61mEE");
00039 amplitudeThreshEB_ = ps.getParameter<double>("amplitudeThresholdEB");
00040 amplitudeThreshEE_ = ps.getParameter<double>("amplitudeThresholdEE");
00041
00042 doEBtimeCorrection_ = ps.getParameter<bool>("doEBtimeCorrection");
00043 doEEtimeCorrection_ = ps.getParameter<bool>("doEEtimeCorrection");
00044 EBtimeCorrAmplitudeBins_ = ps.getParameter<std::vector<double> >("EBtimeCorrAmplitudeBins");
00045 EBtimeCorrShiftBins_ = ps.getParameter<std::vector<double> >("EBtimeCorrShiftBins");
00046 EEtimeCorrAmplitudeBins_ = ps.getParameter<std::vector<double> >("EEtimeCorrAmplitudeBins");
00047 EEtimeCorrShiftBins_ = ps.getParameter<std::vector<double> >("EEtimeCorrShiftBins");
00048 if(EBtimeCorrAmplitudeBins_.size() != EBtimeCorrShiftBins_.size()) {
00049 doEBtimeCorrection_ = false;
00050 edm::LogError("EcalRecHitError") << "Size of EBtimeCorrAmplitudeBins different from EBtimeCorrShiftBins. Forcing no time corrections for EB. ";
00051 }
00052 if(EEtimeCorrAmplitudeBins_.size() != EEtimeCorrShiftBins_.size()) {
00053 doEEtimeCorrection_ = false;
00054 edm::LogError("EcalRecHitError") << "Size of EEtimeCorrAmplitudeBins different from EEtimeCorrShiftBins. Forcing no time corrections for EE. ";
00055 }
00056
00057
00058 ebSpikeThresh_ = ps.getParameter<double>("ebSpikeThreshold");
00059
00060 ebPulseShape_ = ps.getParameter<std::vector<double> >("ebPulseShape");
00061 eePulseShape_ = ps.getParameter<std::vector<double> >("eePulseShape");
00062
00063 kPoorRecoFlagEB_ = ps.getParameter<bool>("kPoorRecoFlagEB");
00064 kPoorRecoFlagEE_ = ps.getParameter<bool>("kPoorRecoFlagEE");;
00065 chi2ThreshEB_=ps.getParameter<double>("chi2ThreshEB_");
00066 chi2ThreshEE_=ps.getParameter<double>("chi2ThreshEE_");
00067 EBchi2Parameters_ = ps.getParameter<std::vector<double> >("EBchi2Parameters");
00068 EEchi2Parameters_ = ps.getParameter<std::vector<double> >("EEchi2Parameters");
00069 }
00070
00071 void
00072 EcalUncalibRecHitWorkerGlobal::set(const edm::EventSetup& es)
00073 {
00074
00075 es.get<EcalGainRatiosRcd>().get(gains);
00076 es.get<EcalPedestalsRcd>().get(peds);
00077
00078
00079 es.get<EcalWeightXtalGroupsRcd>().get(grps);
00080 es.get<EcalTBWeightsRcd>().get(wgts);
00081
00082
00083
00084
00085 es.get<EcalTimeCalibConstantsRcd>().get(itime);
00086 es.get<EcalTimeOffsetConstantRcd>().get(offtime);
00087 }
00088
00089
00090
00091 template < class C >
00092 int EcalUncalibRecHitWorkerGlobal::isSaturated(const C & dataFrame)
00093 {
00094
00095 int cnt;
00096 for (int j = 0; j < C::MAXSAMPLES - 5; ++j) {
00097 cnt = 0;
00098 for (int i = j; i < (j + 5) && i < C::MAXSAMPLES; ++i) {
00099 if ( dataFrame.sample(i).gainId() == 0 ) ++cnt;
00100 }
00101 if ( cnt == 5 ) return j-1 ;
00102 }
00103 return -1;
00104 }
00105
00106
00107 double EcalUncalibRecHitWorkerGlobal::timeCorrectionEB(float ampliEB){
00108
00109 double theCorrection=0;
00110
00111
00112 int myBin = -1;
00113 for (int bin=0; bin<(int)EBtimeCorrAmplitudeBins_.size(); bin++ ){
00114 if(ampliEB > EBtimeCorrAmplitudeBins_.at(bin)) {
00115 myBin = bin; }
00116 else break;
00117 }
00118
00119 if (myBin == -1)
00120 {
00121 theCorrection = EBtimeCorrShiftBins_.at(0);
00122 }
00123 else if ( myBin == ((int)(EBtimeCorrAmplitudeBins_.size()-1)) )
00124 {
00125 theCorrection = EBtimeCorrShiftBins_.at( myBin );
00126 }
00127 else if ( -1 < myBin && myBin < ((int)EBtimeCorrAmplitudeBins_.size()-1) )
00128 {
00129
00130 theCorrection = ( EBtimeCorrShiftBins_.at(myBin+1) - EBtimeCorrShiftBins_.at(myBin) );
00131 theCorrection *= ( ((double)ampliEB) - EBtimeCorrAmplitudeBins_.at(myBin) ) / ( EBtimeCorrAmplitudeBins_.at(myBin+1) - EBtimeCorrAmplitudeBins_.at(myBin) );
00132 theCorrection += EBtimeCorrShiftBins_.at(myBin);
00133 }
00134 else
00135 {
00136 edm::LogError("EcalRecHitError") << "Assigning time correction impossible. Setting it to 0 ";
00137 theCorrection = 0.;
00138 }
00139
00140
00141 return theCorrection/25.;
00142 }
00143
00144
00145 double EcalUncalibRecHitWorkerGlobal::timeCorrectionEE(float ampliEE){
00146
00147 double theCorrection=0;
00148
00149 int myBin = -1;
00150 for (int bin=0; bin<(int)EEtimeCorrAmplitudeBins_.size(); bin++ ){
00151 if(ampliEE > EEtimeCorrAmplitudeBins_.at(bin)) {
00152 myBin = bin; }
00153 else break;
00154 }
00155
00156 if (myBin == -1)
00157 {
00158 theCorrection = EEtimeCorrShiftBins_.at(0);
00159 }
00160 else if ( myBin == ((int)(EEtimeCorrAmplitudeBins_.size()-1)) )
00161 {
00162 theCorrection = EEtimeCorrShiftBins_.at( myBin );
00163 }
00164 else if ( -1 < myBin && myBin < ((int)EEtimeCorrAmplitudeBins_.size()-1) )
00165 {
00166
00167 theCorrection = ( EEtimeCorrShiftBins_.at(myBin+1) - EEtimeCorrShiftBins_.at(myBin) );
00168 theCorrection *= ( ((double)ampliEE) - EEtimeCorrAmplitudeBins_.at(myBin) ) / ( EEtimeCorrAmplitudeBins_.at(myBin+1) - EEtimeCorrAmplitudeBins_.at(myBin) );
00169 theCorrection += EEtimeCorrShiftBins_.at(myBin);
00170 }
00171 else
00172 {
00173 edm::LogError("EcalRecHitError") << "Assigning time correction impossible. Setting it to 0 ";
00174 theCorrection = 0.;
00175 }
00176
00177
00178 return theCorrection/25.;
00179 }
00180
00181
00182
00183
00184 bool
00185 EcalUncalibRecHitWorkerGlobal::run( const edm::Event & evt,
00186 const EcalDigiCollection::const_iterator & itdg,
00187 EcalUncalibratedRecHitCollection & result )
00188 {
00189 DetId detid(itdg->id());
00190
00191
00192 EcalUncalibratedRecHit uncalibRecHit;
00193
00194
00195 const EcalPedestals::Item * aped = 0;
00196 const EcalMGPAGainRatio * aGain = 0;
00197 const EcalXtalGroupId * gid = 0;
00198 float offsetTime = 0;
00199
00200 if (detid.subdetId()==EcalEndcap) {
00201 unsigned int hashedIndex = EEDetId(detid).hashedIndex();
00202 aped = &peds->endcap(hashedIndex);
00203 aGain = &gains->endcap(hashedIndex);
00204 gid = &grps->endcap(hashedIndex);
00205 offsetTime = offtime->getEEValue();
00206 } else {
00207 unsigned int hashedIndex = EBDetId(detid).hashedIndex();
00208 aped = &peds->barrel(hashedIndex);
00209 aGain = &gains->barrel(hashedIndex);
00210 gid = &grps->barrel(hashedIndex);
00211 offsetTime = offtime->getEBValue();
00212 }
00213
00214 pedVec[0] = aped->mean_x12;
00215 pedVec[1] = aped->mean_x6;
00216 pedVec[2] = aped->mean_x1;
00217 pedRMSVec[0] = aped->rms_x12;
00218 pedRMSVec[1] = aped->rms_x6;
00219 pedRMSVec[2] = aped->rms_x1;
00220 gainRatios[0] = 1.;
00221 gainRatios[1] = aGain->gain12Over6();
00222 gainRatios[2] = aGain->gain6Over1()*aGain->gain12Over6();
00223
00224
00225 EcalTimeCalibConstantMap::const_iterator it = itime->find( detid );
00226 EcalTimeCalibConstant itimeconst = 0;
00227 if( it != itime->end() ) {
00228 itimeconst = (*it);
00229 } else {
00230 edm::LogError("EcalRecHitError") << "No time intercalib const found for xtal "
00231 << detid.rawId()
00232 << "! something wrong with EcalTimeCalibConstants in your DB? ";
00233 }
00234
00235
00236
00237 int leadingSample = -1;
00238 if (detid.subdetId()==EcalEndcap) {
00239 leadingSample = ((EcalDataFrame)(*itdg)).lastUnsaturatedSample();
00240 } else {
00241 leadingSample = ((EcalDataFrame)(*itdg)).lastUnsaturatedSample();
00242 }
00243
00244 if ( leadingSample >= 0 ) {
00245 if ( leadingSample != 4 ) {
00246
00247
00248 float sratio = 1;
00249 if ( detid.subdetId()==EcalBarrel) {
00250 sratio = ebPulseShape_[5] / ebPulseShape_[4];
00251 } else {
00252 sratio = eePulseShape_[5] / eePulseShape_[4];
00253 }
00254 uncalibRecHit = EcalUncalibratedRecHit( (*itdg).id(), 4095*12*sratio, 0, 0, 0);
00255 uncalibRecHit.setRecoFlag( EcalUncalibratedRecHit::kSaturated );
00256 } else {
00257
00258
00259 if (detid.subdetId()==EcalEndcap) {
00260 leadingEdgeMethod_endcap_.setPulseShape( eePulseShape_ );
00261
00262
00263
00264
00265 leadingEdgeMethod_endcap_.setLeadingEdgeSample( leadingSample );
00266 uncalibRecHit = leadingEdgeMethod_endcap_.makeRecHit(*itdg, pedVec, gainRatios, 0, 0);
00267 uncalibRecHit.setRecoFlag( EcalUncalibratedRecHit::kLeadingEdgeRecovered );
00268 leadingEdgeMethod_endcap_.setLeadingEdgeSample( -1 );
00269 } else {
00270 leadingEdgeMethod_barrel_.setPulseShape( ebPulseShape_ );
00271
00272
00273
00274
00275 leadingEdgeMethod_barrel_.setLeadingEdgeSample( leadingSample );
00276 uncalibRecHit = leadingEdgeMethod_barrel_.makeRecHit(*itdg, pedVec, gainRatios, 0, 0);
00277 uncalibRecHit.setRecoFlag( EcalUncalibratedRecHit::kLeadingEdgeRecovered );
00278 leadingEdgeMethod_barrel_.setLeadingEdgeSample( -1 );
00279 }
00280 }
00281
00282 uncalibRecHit.setChi2(0);
00283 uncalibRecHit.setOutOfTimeChi2(0);
00284 } else {
00285
00286 EcalTBWeights::EcalTDCId tdcid(1);
00287 EcalTBWeights::EcalTBWeightMap const & wgtsMap = wgts->getMap();
00288 EcalTBWeights::EcalTBWeightMap::const_iterator wit;
00289 wit = wgtsMap.find( std::make_pair(*gid,tdcid) );
00290 if( wit == wgtsMap.end() ) {
00291 edm::LogError("EcalUncalibRecHitError") << "No weights found for EcalGroupId: "
00292 << gid->id() << " and EcalTDCId: " << tdcid
00293 << "\n skipping digi with id: " << detid.rawId();
00294
00295 return false;
00296 }
00297 const EcalWeightSet& wset = wit->second;
00298
00299 const EcalWeightSet::EcalWeightMatrix& mat1 = wset.getWeightsBeforeGainSwitch();
00300 const EcalWeightSet::EcalWeightMatrix& mat2 = wset.getWeightsAfterGainSwitch();
00301
00302 weights[0] = &mat1;
00303 weights[1] = &mat2;
00304
00305
00306 if (detid.subdetId()==EcalEndcap) {
00307 uncalibRecHit = weightsMethod_endcap_.makeRecHit(*itdg, pedVec, pedRMSVec, gainRatios, weights, testbeamEEShape);
00308 } else {
00309 uncalibRecHit = weightsMethod_barrel_.makeRecHit(*itdg, pedVec, pedRMSVec, gainRatios, weights, testbeamEBShape);
00310 }
00311
00312
00313
00314 float const clockToNsConstant = 25.;
00315 if (detid.subdetId()==EcalEndcap) {
00316 ratioMethod_endcap_.init( *itdg, pedVec, pedRMSVec, gainRatios );
00317 ratioMethod_endcap_.computeTime( EEtimeFitParameters_, EEtimeFitLimits_, EEamplitudeFitParameters_ );
00318 ratioMethod_endcap_.computeAmplitude( EEamplitudeFitParameters_);
00319 EcalUncalibRecHitRatioMethodAlgo<EEDataFrame>::CalculatedRecHit crh = ratioMethod_endcap_.getCalculatedRecHit();
00320 double theTimeCorrectionEE=0;
00321 if(doEEtimeCorrection_) theTimeCorrectionEE = timeCorrectionEE( uncalibRecHit.amplitude() );
00322 uncalibRecHit.setJitter( crh.timeMax - 5 + theTimeCorrectionEE);
00323 uncalibRecHit.setJitterError( std::sqrt(pow(crh.timeError,2) + std::pow(EEtimeConstantTerm_,2)/std::pow(clockToNsConstant,2)) );
00324 uncalibRecHit.setOutOfTimeEnergy( crh.amplitudeMax );
00325
00326 if (uncalibRecHit.amplitude() > pedRMSVec[0] * amplitudeThreshEE_){
00327 float outOfTimeThreshP = outOfTimeThreshG12pEE_;
00328 float outOfTimeThreshM = outOfTimeThreshG12mEE_;
00329
00330
00331
00332 if (uncalibRecHit.amplitude() > 3000.){
00333 for (int iSample = 0; iSample < EEDataFrame::MAXSAMPLES; iSample++) {
00334 int GainId = ((EcalDataFrame)(*itdg)).sample(iSample).gainId();
00335 if (GainId!=1) {
00336 outOfTimeThreshP = outOfTimeThreshG61pEE_;
00337 outOfTimeThreshM = outOfTimeThreshG61mEE_;
00338 break;}
00339 }}
00340 float correctedTime = (crh.timeMax-5) * clockToNsConstant + itimeconst + offsetTime;
00341 float cterm = EEtimeConstantTerm_;
00342 float sigmaped = pedRMSVec[0];
00343 float nterm = EEtimeNconst_*sigmaped/uncalibRecHit.amplitude();
00344 float sigmat = std::sqrt( nterm*nterm + cterm*cterm );
00345 if ( ( correctedTime > sigmat*outOfTimeThreshP ) ||
00346 ( correctedTime < (-1.*sigmat*outOfTimeThreshM) ))
00347 { uncalibRecHit.setRecoFlag( EcalUncalibratedRecHit::kOutOfTime ); }
00348 }
00349
00350 } else {
00351 ratioMethod_barrel_.init( *itdg, pedVec, pedRMSVec, gainRatios );
00352 ratioMethod_barrel_.fixMGPAslew(*itdg);
00353 ratioMethod_barrel_.computeTime( EBtimeFitParameters_, EBtimeFitLimits_, EBamplitudeFitParameters_ );
00354 ratioMethod_barrel_.computeAmplitude( EBamplitudeFitParameters_);
00355 EcalUncalibRecHitRatioMethodAlgo<EBDataFrame>::CalculatedRecHit crh = ratioMethod_barrel_.getCalculatedRecHit();
00356 double theTimeCorrectionEB=0;
00357 if(doEBtimeCorrection_) theTimeCorrectionEB = timeCorrectionEB( uncalibRecHit.amplitude() );
00358
00359
00360
00361
00362
00363
00364
00365 uncalibRecHit.setJitter( crh.timeMax - 5 + theTimeCorrectionEB);
00366
00367 uncalibRecHit.setJitterError( std::sqrt(std::pow(crh.timeError,2) + std::pow(EBtimeConstantTerm_,2)/std::pow(clockToNsConstant,2)) );
00368 uncalibRecHit.setOutOfTimeEnergy( crh.amplitudeMax );
00369
00370 if (uncalibRecHit.amplitude() > pedRMSVec[0] * amplitudeThreshEB_){
00371 float outOfTimeThreshP = outOfTimeThreshG12pEB_;
00372 float outOfTimeThreshM = outOfTimeThreshG12mEB_;
00373
00374
00375
00376 if (uncalibRecHit.amplitude() > 3000.){
00377 for (int iSample = 0; iSample < EBDataFrame::MAXSAMPLES; iSample++) {
00378 int GainId = ((EcalDataFrame)(*itdg)).sample(iSample).gainId();
00379 if (GainId!=1) {
00380 outOfTimeThreshP = outOfTimeThreshG61pEB_;
00381 outOfTimeThreshM = outOfTimeThreshG61mEB_;
00382 break;}
00383 } }
00384 float correctedTime = (crh.timeMax-5) * clockToNsConstant + itimeconst + offsetTime;
00385 float cterm = EBtimeConstantTerm_;
00386 float sigmaped = pedRMSVec[0];
00387 float nterm = EBtimeNconst_*sigmaped/uncalibRecHit.amplitude();
00388 float sigmat = std::sqrt( nterm*nterm + cterm*cterm );
00389 if ( ( correctedTime > sigmat*outOfTimeThreshP ) ||
00390 ( correctedTime < (-1.*sigmat*outOfTimeThreshM) ))
00391 { uncalibRecHit.setRecoFlag( EcalUncalibratedRecHit::kOutOfTime ); }
00392 }
00393 }
00394
00395
00396 if (detid.subdetId()==EcalEndcap) {
00397
00398 double amplitude = uncalibRecHit.amplitude();
00399 double amplitudeOutOfTime = uncalibRecHit.outOfTimeEnergy();
00400 double jitter= uncalibRecHit.jitter();
00401
00402
00403
00404 EcalUncalibRecHitRecChi2Algo<EEDataFrame>chi2expressEE_(
00405 *itdg,
00406 amplitude,
00407 (itimeconst + offsetTime),
00408 amplitudeOutOfTime,
00409 jitter,
00410 pedVec,
00411 pedRMSVec,
00412 gainRatios,
00413 testbeamEEShape,
00414 EEchi2Parameters_
00415 );
00416 double chi2 = chi2expressEE_.chi2();
00417 uncalibRecHit.setChi2(chi2);
00418 double chi2OutOfTime = chi2expressEE_.chi2OutOfTime();
00419 uncalibRecHit.setOutOfTimeChi2(chi2OutOfTime);
00420
00421 if(kPoorRecoFlagEE_)
00422 {
00423 if(chi2>chi2ThreshEE_)uncalibRecHit.setRecoFlag(EcalUncalibratedRecHit::kPoorReco);
00424 }
00425
00426 } else {
00427 double amplitude = uncalibRecHit.amplitude();
00428 double amplitudeOutOfTime = uncalibRecHit.outOfTimeEnergy();
00429 double jitter= uncalibRecHit.jitter();
00430
00431 EcalUncalibRecHitRecChi2Algo<EBDataFrame>chi2expressEB_(
00432 *itdg,
00433 amplitude,
00434 (itimeconst + offsetTime),
00435 amplitudeOutOfTime,
00436 jitter,
00437 pedVec,
00438 pedRMSVec,
00439 gainRatios,
00440 testbeamEBShape,
00441 EBchi2Parameters_
00442 );
00443 double chi2 = chi2expressEB_.chi2();
00444 uncalibRecHit.setChi2(chi2);
00445 double chi2OutOfTime = chi2expressEB_.chi2OutOfTime();
00446 uncalibRecHit.setOutOfTimeChi2(chi2OutOfTime);
00447
00448 if(kPoorRecoFlagEB_)
00449 {
00450 if(chi2>chi2ThreshEB_)uncalibRecHit.setRecoFlag(EcalUncalibratedRecHit::kPoorReco);
00451 }
00452 }
00453 }
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465 if (detid.subdetId()==EcalEndcap) {
00466 result.push_back( uncalibRecHit );
00467 } else {
00468 result.push_back( uncalibRecHit );
00469 }
00470 return true;
00471 }
00472
00473 #include "FWCore/Framework/interface/MakerMacros.h"
00474 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitWorkerFactory.h"
00475 DEFINE_EDM_PLUGIN( EcalUncalibRecHitWorkerFactory, EcalUncalibRecHitWorkerGlobal, "EcalUncalibRecHitWorkerGlobal" );