CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_10_patch1/src/RecoLocalCalo/EcalRecAlgos/src/EcalCleaningAlgo.cc

Go to the documentation of this file.
00001 /* Implementation of class EcalCleaningAlgo
00002    \author Stefano Argiro
00003    \version $Id: EcalCleaningAlgo.cc,v 1.9 2011/05/17 12:07:23 argiro Exp $
00004    \date 20 Dec 2010
00005 */    
00006 
00007 
00008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00009 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00010 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00011 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h" 
00012 
00013 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalCleaningAlgo.h"
00014 
00015 EcalCleaningAlgo::EcalCleaningAlgo(const edm::ParameterSet&  p){
00016   
00017   cThreshold_barrel_   = p.getParameter<double>("cThreshold_barrel");;
00018   cThreshold_endcap_   = p.getParameter<double>("cThreshold_endcap");;  
00019   e4e1_a_barrel_       = p.getParameter<double>("e4e1_a_barrel");
00020   e4e1_b_barrel_       = p.getParameter<double>("e4e1_b_barrel");
00021   e4e1_a_endcap_       = p.getParameter<double>("e4e1_a_endcap");
00022   e4e1_b_endcap_       = p.getParameter<double>("e4e1_b_endcap");
00023   e4e1Treshold_barrel_ = p.getParameter<double>("e4e1Threshold_barrel");
00024   e4e1Treshold_endcap_ = p.getParameter<double>("e4e1Threshold_endcap");
00025 
00026   cThreshold_double_   = p.getParameter<double>("cThreshold_double"); 
00027 
00028   ignoreOutOfTimeThresh_   =p.getParameter<double>("ignoreOutOfTimeThresh");  
00029   tightenCrack_e1_single_  =p.getParameter<double>("tightenCrack_e1_single");
00030   tightenCrack_e4e1_single_=p.getParameter<double>("tightenCrack_e4e1_single");
00031   tightenCrack_e1_double_  =p.getParameter<double>("tightenCrack_e1_double");
00032   tightenCrack_e6e2_double_=p.getParameter<double>("tightenCrack_e6e2_double");
00033   e6e2thresh_=              p.getParameter<double>("e6e2thresh");
00034 
00035 
00036   
00037 }
00038 
00039 
00040 
00041 
00069 EcalRecHit::Flags 
00070 EcalCleaningAlgo::checkTopology(const DetId& id,
00071                                 const EcalRecHitCollection& rhs){
00072 
00073 
00074   float a=0,b=0,e4e1thresh=0,ethresh=0;
00075 
00076 
00077   if( id.subdetId() == EcalBarrel) {
00078     a= e4e1_a_barrel_;
00079     b= e4e1_b_barrel_; 
00080     ethresh=cThreshold_barrel_;
00081    
00082   }
00083   else if( id.subdetId() == EcalEndcap){
00084     a= e4e1_a_endcap_;
00085     b= e4e1_b_endcap_;
00086     ethresh=cThreshold_endcap_;
00087    }
00088   
00089   
00090 
00091   // for energies below threshold, we don't apply e4e1 cut
00092   float energy = recHitE(id,rhs,false);
00093  
00094   if (energy< ethresh) return EcalRecHit::kGood;
00095   if (isNearCrack(id) && energy < ethresh*tightenCrack_e1_single_) 
00096     return EcalRecHit::kGood;
00097 
00098 
00099   float e4e1value = e4e1(id,rhs);
00100   e4e1thresh = a* log10(energy) + b;
00101 
00102   // near cracks the cut is tighter by a factor 
00103   if (isNearCrack(id)) {
00104     e4e1thresh/=tightenCrack_e4e1_single_;
00105   }
00106 
00107   // identify spike
00108   if (e4e1value < e4e1thresh) return EcalRecHit::kWeird; 
00109   
00110 
00111 
00112   // now for double spikes
00113  
00114   // no checking for double spikes in EE
00115   if( id.subdetId() == EcalEndcap) return EcalRecHit::kGood;
00116 
00117   float e6e2value = e6e2(id,rhs);
00118   float e6e2thresh = e6e2thresh_ ;
00119   if (isNearCrack(id) && energy < cThreshold_double_ *tightenCrack_e1_double_ )
00120     return EcalRecHit::kGood;
00121 
00122   if  (energy <  cThreshold_double_) return EcalRecHit::kGood;
00123   
00124   // near cracks the cut is tighter by a factor 
00125   if (id.subdetId() == EcalBarrel && isNearCrack(id)) 
00126     e6e2thresh/=tightenCrack_e6e2_double_;
00127 
00128   // identify double spike
00129   if (e6e2value < e6e2thresh) return EcalRecHit::kDiWeird; 
00130 
00131   return EcalRecHit::kGood;
00132 
00133 }
00134 
00135 
00136 
00137 
00138 
00139 
00140 float EcalCleaningAlgo::e4e1(const DetId& id, 
00141                              const EcalRecHitCollection& rhs){
00142 
00143  
00144   float s4 = 0;
00145   float e1 = recHitE( id, rhs, false );
00146   
00147   
00148   if ( e1 == 0 ) return 0;
00149   const std::vector<DetId>& neighs =  neighbours(id);
00150   for (size_t i=0; i<neighs.size(); ++i)
00151     // avoid hits out of time when making s4
00152     s4+=recHitE(neighs[i],rhs, true);
00153   
00154   return s4 / e1;
00155   
00156  
00157 }
00158 
00159 
00160 
00161 
00162 float EcalCleaningAlgo::e6e2(const DetId& id, 
00163                              const EcalRecHitCollection& rhs){
00164 
00165     float s4_1 = 0;
00166     float s4_2 = 0;
00167     float e1 = recHitE( id, rhs , false );
00168 
00169 
00170     float maxene=0;
00171     DetId maxid;
00172 
00173     if ( e1 == 0 ) return 0;
00174 
00175     const std::vector<DetId>& neighs =  neighbours(id);
00176 
00177     // find the most energetic neighbour ignoring time info
00178     for (size_t i=0; i<neighs.size(); ++i){
00179       float ene = recHitE(neighs[i],rhs,false);
00180       if (ene>maxene)  {
00181         maxene=ene;
00182         maxid = neighs[i];
00183       }
00184     }
00185 
00186     float e2=maxene;
00187 
00188     s4_1 = e4e1(id,rhs)* e1;
00189     s4_2 = e4e1(maxid,rhs)* e2;
00190 
00191     return (s4_1 + s4_2) / (e1+e2) -1. ;
00192 
00193 }
00194 
00195 
00196 
00197 
00198 float EcalCleaningAlgo::recHitE( const DetId id, 
00199                                  const EcalRecHitCollection &recHits,
00200                                  bool useTimingInfo )
00201 {
00202   if ( id.rawId() == 0 ) return 0;
00203   
00204 
00205   float threshold = e4e1Treshold_barrel_;
00206   if ( id.subdetId() == EcalEndcap) threshold = e4e1Treshold_endcap_; 
00207 
00208   EcalRecHitCollection::const_iterator it = recHits.find( id );
00209   if ( it != recHits.end() ){
00210     float ene= (*it).energy();
00211 
00212     // ignore out of time in EB when making e4e1 if so configured
00213     if (useTimingInfo){
00214 
00215       if (id.subdetId()==EcalBarrel &&
00216           it->checkFlag(EcalRecHit::kOutOfTime) 
00217           && ene>ignoreOutOfTimeThresh_) return 0;
00218     }
00219 
00220     // ignore hits below threshold
00221     if (ene < threshold) return 0;
00222 
00223     // else return the energy of this hit
00224     return ene;
00225   }
00226   return 0;
00227 }
00228 
00230 const std::vector<DetId> EcalCleaningAlgo::neighbours(const DetId& id){
00231  
00232   std::vector<DetId> ret;
00233 
00234   if ( id.subdetId() == EcalBarrel) {
00235 
00236     ret.push_back( EBDetId::offsetBy( id,  1, 0 ));
00237     ret.push_back( EBDetId::offsetBy( id, -1, 0 ));
00238     ret.push_back( EBDetId::offsetBy( id,  0, 1 ));
00239     ret.push_back( EBDetId::offsetBy( id,  0,-1 ));
00240   }
00241   // nobody understands what polymorphism is for, sgrunt !
00242   else  if (id.subdetId() == EcalEndcap) {
00243     ret.push_back( EEDetId::offsetBy( id,  1, 0 ));
00244     ret.push_back( EEDetId::offsetBy( id, -1, 0 ));
00245     ret.push_back( EEDetId::offsetBy( id,  0, 1 ));
00246     ret.push_back( EEDetId::offsetBy( id,  0,-1 ));
00247 
00248   }
00249 
00250 
00251   return ret;
00252 
00253 } 
00254 
00255 
00256 bool EcalCleaningAlgo::isNearCrack(const DetId& id){
00257 
00258   if (id.subdetId() == EcalEndcap) { 
00259     return EEDetId::isNextToRingBoundary(id);
00260   } else {
00261     return EBDetId::isNextToBoundary(id);
00262   }
00263 }
00264 
00265 
00266 
00267 void EcalCleaningAlgo::setFlags(EcalRecHitCollection& rhs){
00268   EcalRecHitCollection::iterator rh;
00269   //changing the collection on place
00270   for (rh=rhs.begin(); rh!=rhs.end(); ++rh){
00271     EcalRecHit::Flags state=checkTopology(rh->id(),rhs);
00272     if (state!=EcalRecHit::kGood) { 
00273       rh->unsetFlag(EcalRecHit::kGood);
00274       rh->setFlag(state);
00275     }
00276   }
00277 }
00278