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