CMS 3D CMS Logo

CMSSW_4_4_3_patch1/src/RecoLocalCalo/EcalRecProducers/plugins/EcalUncalibRecHitWorkerGlobal.cc

Go to the documentation of this file.
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         // ratio method parameters
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         // amplitude-dependent correction of time
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         // spike threshold
00058         ebSpikeThresh_ = ps.getParameter<double>("ebSpikeThreshold");
00059         // leading edge parameters
00060         ebPulseShape_ = ps.getParameter<std::vector<double> >("ebPulseShape");
00061         eePulseShape_ = ps.getParameter<std::vector<double> >("eePulseShape");
00062         // chi2 parameters
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         // common setup
00075         es.get<EcalGainRatiosRcd>().get(gains);
00076         es.get<EcalPedestalsRcd>().get(peds);
00077 
00078         // for the weights method
00079         es.get<EcalWeightXtalGroupsRcd>().get(grps);
00080         es.get<EcalTBWeightsRcd>().get(wgts);
00081 
00082         // for the ratio method
00083 
00084         // for the leading edge method
00085         es.get<EcalTimeCalibConstantsRcd>().get(itime);
00086         es.get<EcalTimeOffsetConstantRcd>().get(offtime);
00087 }
00088 
00089 
00090 // check saturation: 5 samples with gainId = 0
00091 template < class C >
00092 int EcalUncalibRecHitWorkerGlobal::isSaturated(const C & dataFrame)
00093 {
00094         //bool saturated_ = 0;
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 ; // the last unsaturated sample
00102         }
00103         return -1; // no saturation found
00104 }
00105 
00106 
00107 double EcalUncalibRecHitWorkerGlobal::timeCorrectionEB(float ampliEB){
00108   // computed initially in ns. Than turned in the BX's, as EcalUncalibratedRecHit need be.
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       // interpolate linearly between two assingned points
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   // convert ns into clocks
00141   return theCorrection/25.;
00142 }
00143 
00144 
00145 double EcalUncalibRecHitWorkerGlobal::timeCorrectionEE(float ampliEE){
00146   // computed initially in ns. Than turned in the BX's, as EcalUncalibratedRecHit need be.
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       // interpolate linearly between two assingned points
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   // convert ns into clocks
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         // intelligence for recHit computation
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         // compute the right bin of the pulse shape using time calibration constants
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         // === amplitude computation ===
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 ) { // saturation
00245                 if ( leadingSample != 4 ) {
00246                         // all samples different from the fifth are not reliable for the amplitude estimation
00247                         // put by default the energy at the saturation threshold and flag as saturated
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                         // float clockToNsConstant = 25.;
00258                         // reconstruct the rechit
00259                         if (detid.subdetId()==EcalEndcap) {
00260                                 leadingEdgeMethod_endcap_.setPulseShape( eePulseShape_ );
00261                                 // float mult = (float)eePulseShape_.size() / (float)(*itdg).size();
00262                                 // bin (or some analogous mapping) will be used instead of the leadingSample
00263                                 //int bin  = (int)(( (mult * leadingSample + mult/2) * clockToNsConstant + itimeconst ) / clockToNsConstant);
00264                                 // bin is not uset for the moment
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                                 // float mult = (float)ebPulseShape_.size() / (float)(*itdg).size();
00272                                 // bin (or some analogous mapping) will be used instead of the leadingSample
00273                                 //int bin  = (int)(( (mult * leadingSample + mult/2) * clockToNsConstant + itimeconst ) / clockToNsConstant);
00274                                 // bin is not uset for the moment
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);         // do not propagate the default chi2 = -1 value to the calib rechit (mapped to 64), set it to 0 when saturation
00283                 uncalibRecHit.setOutOfTimeChi2(0);
00284         } else {
00285                 // weights method
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; // this is the EcalWeightSet
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                 // get uncalibrated recHit from weights
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                 // === time computation ===
00313                 // ratio method
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                                 // consider flagging as kOutOfTime only if above noise
00326                                 if (uncalibRecHit.amplitude() > pedRMSVec[0] * amplitudeThreshEE_){
00327                                   float outOfTimeThreshP = outOfTimeThreshG12pEE_;
00328                                   float outOfTimeThreshM = outOfTimeThreshG12mEE_;
00329                                   // determine if gain has switched away from gainId==1 (x12 gain)
00330                                   // and determine cuts (number of 'sigmas') to ose for kOutOfTime
00331                                   // >3k ADC is necessasry condition for gain switch to occur
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];  // approx for lower gains
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                                 // the correction for gain switch (when the sample before is ignored in ratioAlgo) is now included in the configurable correction
00359                                 //      bool gainSwitch = ratioMethod_barrel_.fixMGPAslew(*itdg);
00360                                 //      if(gainSwitch){
00361                                 //        uncalibRecHit.setJitter( crh.timeMax - 5 - 0.04 + theTimeCorrectionEB);  // introduce additional 1ns shift
00362                                 //      }else{
00363                                 //        uncalibRecHit.setJitter( crh.timeMax - 5 + theTimeCorrectionEB);
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                                 // consider flagging as kOutOfTime only if above noise
00370                                 if (uncalibRecHit.amplitude() > pedRMSVec[0] * amplitudeThreshEB_){
00371                                   float outOfTimeThreshP = outOfTimeThreshG12pEB_;
00372                                   float outOfTimeThreshM = outOfTimeThreshG12mEB_;
00373                                   // determine if gain has switched away from gainId==1 (x12 gain)
00374                                   // and determine cuts (number of 'sigmas') to ose for kOutOfTime
00375                                   // >3k ADC is necessasry condition for gain switch to occur
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; //GF add here
00385                                   float cterm         = EBtimeConstantTerm_;
00386                                   float sigmaped      = pedRMSVec[0];  // approx for lower gains
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                 // === chi2express ===
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         // remove setting of kFake, which can be misleading for the time being
00456         //if ( detid.subdetId()==EcalBarrel ) {
00457         //        if ( uncalibRecHit.jitter()*25. > -5 ) {
00458         //                EBDataFrame dt(*itdg);
00459         //                if ( dt.spikeEstimator() < ebSpikeThresh_ ) uncalibRecHit.setRecoFlag( EcalUncalibratedRecHit::kFake );
00460         //        }
00461         //}
00462         
00463 
00464         // put the recHit in the collection
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" );