EcalCleaningAlgo.cc

Go to the documentation of this file.

/* Implementation of class EcalCleaningAlgo
   \author Stefano Argiro
   \date 20 Dec 2010
*/    


#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHit.h" 

#include "RecoLocalCalo/EcalRecAlgos/interface/EcalCleaningAlgo.h"

EcalCleaningAlgo::EcalCleaningAlgo(const edm::ParameterSet&  p){
  
  cThreshold_barrel_   = p.getParameter<double>("cThreshold_barrel");;
  cThreshold_endcap_   = p.getParameter<double>("cThreshold_endcap");;  
  e4e1_a_barrel_       = p.getParameter<double>("e4e1_a_barrel");
  e4e1_b_barrel_       = p.getParameter<double>("e4e1_b_barrel");
  e4e1_a_endcap_       = p.getParameter<double>("e4e1_a_endcap");
  e4e1_b_endcap_       = p.getParameter<double>("e4e1_b_endcap");
  e4e1Treshold_barrel_ = p.getParameter<double>("e4e1Threshold_barrel");
  e4e1Treshold_endcap_ = p.getParameter<double>("e4e1Threshold_endcap");

  cThreshold_double_   = p.getParameter<double>("cThreshold_double"); 

  ignoreOutOfTimeThresh_   =p.getParameter<double>("ignoreOutOfTimeThresh");  
  tightenCrack_e1_single_  =p.getParameter<double>("tightenCrack_e1_single");
  tightenCrack_e4e1_single_=p.getParameter<double>("tightenCrack_e4e1_single");
  tightenCrack_e1_double_  =p.getParameter<double>("tightenCrack_e1_double");
  tightenCrack_e6e2_double_=p.getParameter<double>("tightenCrack_e6e2_double");
  e6e2thresh_=              p.getParameter<double>("e6e2thresh");


  
}




EcalRecHit::Flags 
EcalCleaningAlgo::checkTopology(const DetId& id,
                const EcalRecHitCollection& rhs){


  float a=0,b=0,e4e1thresh=0,ethresh=0;


  if( id.subdetId() == EcalBarrel) {
    a= e4e1_a_barrel_;
    b= e4e1_b_barrel_; 
    ethresh=cThreshold_barrel_;
   
  }
  else if( id.subdetId() == EcalEndcap){
    a= e4e1_a_endcap_;
    b= e4e1_b_endcap_;
    ethresh=cThreshold_endcap_;
   }
  
  

  // for energies below threshold, we don't apply e4e1 cut
  float energy = recHitE(id,rhs,false);
 
  if (energy< ethresh) return EcalRecHit::kGood;
  if (isNearCrack(id) && energy < ethresh*tightenCrack_e1_single_) 
    return EcalRecHit::kGood;


  float e4e1value = e4e1(id,rhs);
  e4e1thresh = a* log10(energy) + b;

  // near cracks the cut is tighter by a factor 
  if (isNearCrack(id)) {
    e4e1thresh/=tightenCrack_e4e1_single_;
  }

  // identify spike
  if (e4e1value < e4e1thresh) return EcalRecHit::kWeird; 
  


  // now for double spikes
 
  // no checking for double spikes in EE
  if( id.subdetId() == EcalEndcap) return EcalRecHit::kGood;

  float e6e2value = e6e2(id,rhs);
  float e6e2thresh = e6e2thresh_ ;
  if (isNearCrack(id) && energy < cThreshold_double_ *tightenCrack_e1_double_ )
    return EcalRecHit::kGood;

  if  (energy <  cThreshold_double_) return EcalRecHit::kGood;
  
  // near cracks the cut is tighter by a factor 
  if (id.subdetId() == EcalBarrel && isNearCrack(id)) 
    e6e2thresh/=tightenCrack_e6e2_double_;

  // identify double spike
  if (e6e2value < e6e2thresh) return EcalRecHit::kDiWeird; 

  return EcalRecHit::kGood;

}






float EcalCleaningAlgo::e4e1(const DetId& id, 
                 const EcalRecHitCollection& rhs){

 
  float s4 = 0;
  float e1 = recHitE( id, rhs, false );
  
  
  if ( e1 == 0 ) return 0;
  const std::vector<DetId>& neighs =  neighbours(id);
  for (size_t i=0; i<neighs.size(); ++i)
    // avoid hits out of time when making s4
    s4+=recHitE(neighs[i],rhs, true);
  
  return s4 / e1;
  
 
}




float EcalCleaningAlgo::e6e2(const DetId& id, 
                 const EcalRecHitCollection& rhs){

    float s4_1 = 0;
    float s4_2 = 0;
    float e1 = recHitE( id, rhs , false );


    float maxene=0;
    DetId maxid;

    if ( e1 == 0 ) return 0;

    const std::vector<DetId>& neighs =  neighbours(id);

    // find the most energetic neighbour ignoring time info
    for (size_t i=0; i<neighs.size(); ++i){
      float ene = recHitE(neighs[i],rhs,false);
      if (ene>maxene)  {
    maxene=ene;
    maxid = neighs[i];
      }
    }

    float e2=maxene;

    s4_1 = e4e1(id,rhs)* e1;
    s4_2 = e4e1(maxid,rhs)* e2;

    return (s4_1 + s4_2) / (e1+e2) -1. ;

}




float EcalCleaningAlgo::recHitE( const DetId id, 
                 const EcalRecHitCollection &recHits,
                                 bool useTimingInfo )
{
  if ( id.rawId() == 0 ) return 0;
  

  float threshold = e4e1Treshold_barrel_;
  if ( id.subdetId() == EcalEndcap) threshold = e4e1Treshold_endcap_; 

  EcalRecHitCollection::const_iterator it = recHits.find( id );
  if ( it != recHits.end() ){
    float ene= (*it).energy();

    // ignore out of time in EB when making e4e1 if so configured
    if (useTimingInfo){

      if (id.subdetId()==EcalBarrel &&
      it->checkFlag(EcalRecHit::kOutOfTime) 
      && ene>ignoreOutOfTimeThresh_) return 0;
    }

    // ignore hits below threshold
    if (ene < threshold) return 0;

    // else return the energy of this hit
    return ene;
  }
  return 0;
}

const std::vector<DetId> EcalCleaningAlgo::neighbours(const DetId& id){
 
  std::vector<DetId> ret;

  if ( id.subdetId() == EcalBarrel) {

    ret.push_back( EBDetId::offsetBy( id,  1, 0 ));
    ret.push_back( EBDetId::offsetBy( id, -1, 0 ));
    ret.push_back( EBDetId::offsetBy( id,  0, 1 ));
    ret.push_back( EBDetId::offsetBy( id,  0,-1 ));
  }
  // nobody understands what polymorphism is for, sgrunt !
  else  if (id.subdetId() == EcalEndcap) {
    ret.push_back( EEDetId::offsetBy( id,  1, 0 ));
    ret.push_back( EEDetId::offsetBy( id, -1, 0 ));
    ret.push_back( EEDetId::offsetBy( id,  0, 1 ));
    ret.push_back( EEDetId::offsetBy( id,  0,-1 ));

  }


  return ret;

} 


bool EcalCleaningAlgo::isNearCrack(const DetId& id){

  if (id.subdetId() == EcalEndcap) { 
    return EEDetId::isNextToRingBoundary(id);
  } else {
    return EBDetId::isNextToBoundary(id);
  }
}



void EcalCleaningAlgo::setFlags(EcalRecHitCollection& rhs){
  EcalRecHitCollection::iterator rh;
  //changing the collection on place
  for (rh=rhs.begin(); rh!=rhs.end(); ++rh){
    EcalRecHit::Flags state=checkTopology(rh->id(),rhs);
    if (state!=EcalRecHit::kGood) { 
      rh->unsetFlag(EcalRecHit::kGood);
      rh->setFlag(state);
    }
  }
}