Producer for particle flow rechits (PFRecHit) in HCAL. More...

#include <PFRecHitProducerHCAL.h>

Inheritance diagram for PFRecHitProducerHCAL:

Public Member Functions
	PFRecHitProducerHCAL (const edm::ParameterSet &)

	~PFRecHitProducerHCAL ()

Public Member Functions inherited from PFRecHitProducer
virtual void	beginRun (edm::Run &run, const edm::EventSetup &es)

virtual void	endRun ()

	PFRecHitProducer (const edm::ParameterSet &)

void	produce (edm::Event &iEvent, const edm::EventSetup &iSetup)

	~PFRecHitProducer ()

Public Member Functions inherited from edm::EDProducer
	EDProducer ()

virtual	~EDProducer ()

Public Member Functions inherited from edm::ProducerBase
	ProducerBase ()

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

boost::function< void(const BranchDescription &)>	registrationCallback () const
	used by the fwk to register list of products More...

virtual	~ProducerBase ()

Private Member Functions
reco::PFRecHit *	createHcalRecHit (const DetId &detid, double energy, PFLayer::Layer layer, const CaloSubdetectorGeometry *geom, unsigned newDetId=0)

void	createRecHits (std::vector< reco::PFRecHit > &rechits, std::vector< reco::PFRecHit > &rechitsCleaned, edm::Event &, const edm::EventSetup &)

void	findRecHitNeighbours (reco::PFRecHit &rh, const std::map< unsigned, unsigned > &sortedHits, const CaloSubdetectorTopology &barrelTopo, const CaloSubdetectorGeometry &barrelGeom, const CaloSubdetectorTopology &endcapTopo, const CaloSubdetectorGeometry &endcapGeom)

void	findRecHitNeighboursCT (reco::PFRecHit &rh, const std::map< unsigned, unsigned > &sortedHits, const CaloSubdetectorTopology &topology)

DetId	getNorth (const DetId &id, const CaloSubdetectorTopology &topology)

DetId	getSouth (const DetId &id, const CaloSubdetectorTopology &topology)

Private Attributes
bool	applyLongShortDPG_

bool	applyPulseDPG_

bool	applyTimeDPG_

bool	ECAL_Compensate_

double	ECAL_Compensation_

unsigned int	ECAL_Dead_Code_

double	ECAL_Threshold_

double	EM_Depth_

double	HAD_Depth_

bool	HCAL_Calib_

float	HCAL_Calib_29

int	hcalHFDigiTimeFlagValue_

int	hcalHFInTimeWindowFlagValue_

int	hcalHFLongShortFlagValue_

int	HcalMaxAllowedChannelStatusSev_

int	HcalMaxAllowedHFDigiTimeSev_

int	HcalMaxAllowedHFInTimeWindowSev_

int	HcalMaxAllowedHFLongShortSev_

bool	HF_Calib_

float	HF_Calib_29

edm::InputTag	inputTagCaloTowers_

edm::InputTag	inputTagHcalRecHitsHBHE_

edm::InputTag	inputTagHcalRecHitsHF_

double	longFibre_Cut

double	longFibre_Fraction

double	longShortFibre_Cut

double	maxLongTiming_Cut

double	maxShortTiming_Cut

double	minLongTiming_Cut

double	minShortTiming_Cut

bool	navigation_HF_

double	shortFibre_Cut

double	shortFibre_Fraction

double	thresh_HF_
	threshold for HF More...

double	weight_HFem_

double	weight_HFhad_

Additional Inherited Members
Public Types inherited from edm::EDProducer
typedef EDProducer	ModuleType

typedef WorkerT< EDProducer >	WorkerType

Public Types inherited from edm::ProducerBase
typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Static Public Member Functions inherited from edm::EDProducer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Protected Types inherited from PFRecHitProducer
typedef std::map< unsigned, unsigned >::const_iterator	IDH

Protected Member Functions inherited from edm::EDProducer
CurrentProcessingContext const *	currentContext () const

Protected Member Functions inherited from edm::ProducerBase
template<class TProducer , class TMethod >
void	callWhenNewProductsRegistered (TProducer *iProd, TMethod iMethod)

Protected Attributes inherited from PFRecHitProducer
const HcalPFCorrs *	myPFCorr

const EcalChannelStatus *	theEcalChStatus

const HcalChannelQuality *	theHcalChStatus

const CaloTowerConstituentsMap *	theTowerConstituentsMap

double	thresh_Barrel_
	rechits with E < threshold will not give rise to a PFRecHit More...

double	thresh_Endcap_

bool	verbose_
	verbose ? More...

Detailed Description

Producer for particle flow rechits (PFRecHit) in HCAL.

Author: Colin Bernet

Date: february 2008

Definition at line 33 of file PFRecHitProducerHCAL.h.

Constructor & Destructor Documentation

PFRecHitProducerHCAL::PFRecHitProducerHCAL ( const edm::ParameterSet & iConfig )

explicit

Definition at line 42 of file PFRecHitProducerHCAL.cc.

References applyLongShortDPG_, applyPulseDPG_, applyTimeDPG_, ECAL_Compensate_, ECAL_Compensation_, ECAL_Dead_Code_, ECAL_Threshold_, EM_Depth_, edm::ParameterSet::getParameter(), HAD_Depth_, HCAL_Calib_, HCAL_Calib_29, hcalHFDigiTimeFlagValue_, hcalHFInTimeWindowFlagValue_, hcalHFLongShortFlagValue_, HcalMaxAllowedChannelStatusSev_, HcalMaxAllowedHFDigiTimeSev_, HcalMaxAllowedHFInTimeWindowSev_, HcalMaxAllowedHFLongShortSev_, HF_Calib_, HF_Calib_29, HcalCaloFlagLabels::HFDigiTime, HcalCaloFlagLabels::HFInTimeWindow, HcalCaloFlagLabels::HFLongShort, inputTagCaloTowers_, inputTagHcalRecHitsHBHE_, inputTagHcalRecHitsHF_, longFibre_Cut, longFibre_Fraction, longShortFibre_Cut, maxLongTiming_Cut, maxShortTiming_Cut, minLongTiming_Cut, minShortTiming_Cut, navigation_HF_, shortFibre_Cut, shortFibre_Fraction, thresh_HF_, weight_HFem_, and weight_HFhad_.

   : PFRecHitProducer( iConfig ) 
 {
 
  
 
   // access to the collections of rechits 
 
   
   inputTagHcalRecHitsHBHE_ =
     iConfig.getParameter<InputTag>("hcalRecHitsHBHE");
     
   inputTagHcalRecHitsHF_ =
     iConfig.getParameter<InputTag>("hcalRecHitsHF");
     
  
   inputTagCaloTowers_ = 
     iConfig.getParameter<InputTag>("caloTowers");
    
   thresh_HF_ = 
     iConfig.getParameter<double>("thresh_HF");
   navigation_HF_ = 
     iConfig.getParameter<bool>("navigation_HF");
   weight_HFem_ =
     iConfig.getParameter<double>("weight_HFem");
   weight_HFhad_ =
     iConfig.getParameter<double>("weight_HFhad");
 
   HCAL_Calib_ =
     iConfig.getParameter<bool>("HCAL_Calib");
   HF_Calib_ =
     iConfig.getParameter<bool>("HF_Calib");
   HCAL_Calib_29 = 
     iConfig.getParameter<double>("HCAL_Calib_29");
   HF_Calib_29 = 
     iConfig.getParameter<double>("HF_Calib_29");
 
   shortFibre_Cut = iConfig.getParameter<double>("ShortFibre_Cut");
   longFibre_Fraction = iConfig.getParameter<double>("LongFibre_Fraction");
 
   longFibre_Cut = iConfig.getParameter<double>("LongFibre_Cut");
   shortFibre_Fraction = iConfig.getParameter<double>("ShortFibre_Fraction");
 
   applyLongShortDPG_ = iConfig.getParameter<bool>("ApplyLongShortDPG");
 
   longShortFibre_Cut = iConfig.getParameter<double>("LongShortFibre_Cut");
   minShortTiming_Cut = iConfig.getParameter<double>("MinShortTiming_Cut");
   maxShortTiming_Cut = iConfig.getParameter<double>("MaxShortTiming_Cut");
   minLongTiming_Cut = iConfig.getParameter<double>("MinLongTiming_Cut");
   maxLongTiming_Cut = iConfig.getParameter<double>("MaxLongTiming_Cut");
 
   applyTimeDPG_ = iConfig.getParameter<bool>("ApplyTimeDPG");
   applyPulseDPG_ = iConfig.getParameter<bool>("ApplyPulseDPG");
   HcalMaxAllowedHFLongShortSev_ = iConfig.getParameter<int>("HcalMaxAllowedHFLongShortSev");
   HcalMaxAllowedHFDigiTimeSev_ = iConfig.getParameter<int>("HcalMaxAllowedHFDigiTimeSev");
   HcalMaxAllowedHFInTimeWindowSev_ = iConfig.getParameter<int>("HcalMaxAllowedHFInTimeWindowSev");
   HcalMaxAllowedChannelStatusSev_ = iConfig.getParameter<int>("HcalMaxAllowedChannelStatusSev");
 
   ECAL_Compensate_ = iConfig.getParameter<bool>("ECAL_Compensate");
   ECAL_Threshold_ = iConfig.getParameter<double>("ECAL_Threshold");
   ECAL_Compensation_ = iConfig.getParameter<double>("ECAL_Compensation");
   ECAL_Dead_Code_ = iConfig.getParameter<unsigned int>("ECAL_Dead_Code");
 
   EM_Depth_ = iConfig.getParameter<double>("EM_Depth");
   HAD_Depth_ = iConfig.getParameter<double>("HAD_Depth");
 
   //Get integer values of individual HCAL HF flags
   hcalHFLongShortFlagValue_=1<<HcalCaloFlagLabels::HFLongShort;
   hcalHFDigiTimeFlagValue_=1<<HcalCaloFlagLabels::HFDigiTime;
   hcalHFInTimeWindowFlagValue_=1<<HcalCaloFlagLabels::HFInTimeWindow;
 
 
   //--ab
   produces<reco::PFRecHitCollection>("HFHAD").setBranchAlias("HFHADRecHits");
   produces<reco::PFRecHitCollection>("HFEM").setBranchAlias("HFEMRecHits");
   //--ab
 }

PFRecHitProducerHCAL::~PFRecHitProducerHCAL ( )

Definition at line 122 of file PFRecHitProducerHCAL.cc.

122 {}

Member Function Documentation

reco::PFRecHit * PFRecHitProducerHCAL::createHcalRecHit	(	const DetId &	detid,
		double	energy,
		PFLayer::Layer	layer,
		const CaloSubdetectorGeometry *	geom,
		unsigned	newDetId = `0`
	)

private

Definition at line 1061 of file PFRecHitProducerHCAL.cc.

References cond::rpcobgas::detid, EM_Depth_, CaloCellGeometry::getCorners(), CaloSubdetectorGeometry::getGeometry(), CaloCellGeometry::getPosition(), HAD_Depth_, PFLayer::HF_EM, PFLayer::HF_HAD, position, DetId::rawId(), reco::PFRecHit::setNECorner(), reco::PFRecHit::setNWCorner(), reco::PFRecHit::setSECorner(), reco::PFRecHit::setSWCorner(), x, PV3DBase< T, PVType, FrameType >::x(), detailsBasic3DVector::y, PV3DBase< T, PVType, FrameType >::y(), detailsBasic3DVector::z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by createRecHits().

                                         {
   
   const CaloCellGeometry *thisCell = geom->getGeometry(detid);
   if(!thisCell) {
     edm::LogError("PFRecHitProducerHCAL")
       <<"warning detid "<<detid.rawId()<<" not found in layer "
       <<layer<<endl;
     return 0;
   }
   
   const GlobalPoint& position = thisCell->getPosition();
   
   double depth_correction = 0.;
   switch ( layer ) { 
   case PFLayer::HF_EM:
     depth_correction = position.z() > 0. ? EM_Depth_ : -EM_Depth_;
     break;
   case PFLayer::HF_HAD:
     depth_correction = position.z() > 0. ? HAD_Depth_ : -HAD_Depth_;
     break;
   default:
     break;
   }
 
   unsigned id = detid;
   if(newDetId) id = newDetId;
   reco::PFRecHit *rh = 
     new reco::PFRecHit( id,  layer, energy, 
             position.x(), position.y(), position.z()+depth_correction, 
             0,0,0 );
  
   
   
   
   // set the corners
   const CaloCellGeometry::CornersVec& corners = thisCell->getCorners();
 
   assert( corners.size() == 8 );
 
   rh->setNECorner( corners[0].x(), corners[0].y(),  corners[0].z()+depth_correction );
   rh->setSECorner( corners[1].x(), corners[1].y(),  corners[1].z()+depth_correction );
   rh->setSWCorner( corners[2].x(), corners[2].y(),  corners[2].z()+depth_correction );
   rh->setNWCorner( corners[3].x(), corners[3].y(),  corners[3].z()+depth_correction );
  
   return rh;
 }

void PFRecHitProducerHCAL::createRecHits	(	std::vector< reco::PFRecHit > &	rechits,
		std::vector< reco::PFRecHit > &	rechitsCleaned,
		edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

privatevirtual

gets hcal barrel and endcap rechits, translate them to PFRecHits, which are stored in the rechits vector

if ( !theStatusValue )

Implements PFRecHitProducer.

Definition at line 126 of file PFRecHitProducerHCAL.cc.

References abs, applyLongShortDPG_, applyPulseDPG_, applyTimeDPG_, HiRecoJets_cff::caloTowers, CaloTower::constituents(), CaloTowerConstituentsMap::constituentsOf(), CaloTower::constituentsSize(), createHcalRecHit(), CaloRecHit::detid(), cond::rpcobgas::detid, alignCSCRings::e, DetId::Ecal, ECAL_Compensate_, ECAL_Compensation_, ECAL_Dead_Code_, ECAL_Threshold_, CaloTower::emEnergy(), EcalCondObjectContainer< T >::end(), CaloRecHit::energy(), relval_parameters_module::energy, edm::hlt::Exception, EcalCondObjectContainer< T >::find(), findRecHitNeighbours(), findRecHitNeighboursCT(), newFWLiteAna::found, edm::EventSetup::get(), edm::Event::getByLabel(), HcalSeverityLevelComputer::getSeverityLevel(), HcalChannelStatus::getValue(), HcalCondObjectContainer< Item >::getValues(), CaloTower::hadEnergy(), patZpeak::handle, DetId::Hcal, PFLayer::HCAL_BARREL1, HCAL_Calib_, HCAL_Calib_29, PFLayer::HCAL_ENDCAP, HcalBarrel, HcalEndcap, HcalForward, hcalHFDigiTimeFlagValue_, hcalHFInTimeWindowFlagValue_, hcalHFLongShortFlagValue_, HcalMaxAllowedChannelStatusSev_, HcalMaxAllowedHFDigiTimeSev_, HcalMaxAllowedHFInTimeWindowSev_, HcalMaxAllowedHFLongShortSev_, HcalOuter, HF_Calib_, HF_Calib_29, PFLayer::HF_EM, PFLayer::HF_HAD, i, CaloTower::id(), HcalDetId::ieta(), inputTagCaloTowers_, inputTagHcalRecHitsHBHE_, inputTagHcalRecHitsHF_, HcalDetId::iphi(), edm::HandleBase::isValid(), j, longFibre_Cut, longFibre_Fraction, longShortFibre_Cut, max(), maxLongTiming_Cut, maxShortTiming_Cut, minLongTiming_Cut, minShortTiming_Cut, edm::ESHandle< class >::product(), edm::Event::put(), DetId::rawId(), reco::PFRecHit::setEnergyUp(), reco::PFRecHit::setRescale(), shortFibre_Cut, shortFibre_Fraction, HcalDetId::subdet(), PFRecHitProducer::theEcalChStatus, PFRecHitProducer::theHcalChStatus, PFRecHitProducer::theTowerConstituentsMap, PFRecHitProducer::thresh_Barrel_, PFRecHitProducer::thresh_Endcap_, thresh_HF_, weight_HFem_, and weight_HFhad_.

                                                    {
 
   
   // this map is necessary to find the rechit neighbours efficiently
   //C but I should think about using Florian's hashed index to do this.
   //C in which case the map might not be necessary anymore
   //C however the hashed index does not seem to be implemented for HCAL
   // 
   // the key of this map is detId. 
   // the value is the index in the rechits vector
   map<unsigned,  unsigned > idSortedRecHits;
   map<unsigned,  unsigned > idSortedRecHitsHFEM;
   map<unsigned,  unsigned > idSortedRecHitsHFHAD;
   typedef map<unsigned, unsigned >::iterator IDH;  
 
 
   edm::ESHandle<CaloGeometry> geoHandle;
   iSetup.get<CaloGeometryRecord>().get(geoHandle);
   
   // get the hcalBarrel geometry
   const CaloSubdetectorGeometry *hcalBarrelGeometry = 
     geoHandle->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
 
   // get the endcap geometry
   const CaloSubdetectorGeometry *hcalEndcapGeometry = 
     geoHandle->getSubdetectorGeometry(DetId::Hcal, HcalEndcap);
 
   // Get Hcal Severity Level Computer, so that the severity of each rechit flag/status may be determined
   edm::ESHandle<HcalSeverityLevelComputer> hcalSevLvlComputerHndl;
   iSetup.get<HcalSeverityLevelComputerRcd>().get(hcalSevLvlComputerHndl);
   const HcalSeverityLevelComputer* hcalSevLvlComputer = hcalSevLvlComputerHndl.product();
 
   //--ab
   auto_ptr< vector<reco::PFRecHit> > HFHADRecHits( new vector<reco::PFRecHit> ); 
   auto_ptr< vector<reco::PFRecHit> > HFEMRecHits( new vector<reco::PFRecHit> ); 
   //--ab
 
   // 2 possibilities to make HCAL clustering :
   // - from the HCAL rechits
   // - from the CaloTowers. 
   // ultimately, clustering will be done taking CaloTowers as an 
   // input. This possibility is currently under investigation, and 
   // was thus made optional.
 
   // in the first step, we will fill the map of PFRecHits hcalrechits
   // either from CaloTowers or from HCAL rechits. 
 
   // in the second step, we will perform clustering on this map.
 
   if( !(inputTagCaloTowers_ == InputTag()) ) {
       
     edm::Handle<CaloTowerCollection> caloTowers; 
     CaloTowerTopology caloTowerTopology; 
     const CaloSubdetectorGeometry *caloTowerGeometry = 0; 
     // = geometry_->getSubdetectorGeometry(id)
 
     // get calotowers
     bool found = iEvent.getByLabel(inputTagCaloTowers_,
                    caloTowers);
 
     if(!found) {
       ostringstream err;
       err<<"could not find rechits "<<inputTagCaloTowers_;
       LogError("PFRecHitProducerHCAL")<<err.str()<<endl;
     
       throw cms::Exception( "MissingProduct", err.str());
     }
     else {
       assert( caloTowers.isValid() );
 
       // get HF rechits
       edm::Handle<HFRecHitCollection>  hfHandle;  
       found = iEvent.getByLabel(inputTagHcalRecHitsHF_,
                 hfHandle);
       
       if(!found) {
     ostringstream err;
     err<<"could not find HF rechits "<<inputTagHcalRecHitsHF_;
     LogError("PFRecHitProducerHCAL")<<err.str()<<endl;
     
     throw cms::Exception( "MissingProduct", err.str());
       }
       else {
     assert( hfHandle.isValid() );
       }
       
       // get HBHE rechits
       edm::Handle<HBHERecHitCollection>  hbheHandle;  
       found = iEvent.getByLabel(inputTagHcalRecHitsHBHE_,
                 hbheHandle);
       
       if(!found) {
     ostringstream err;
     err<<"could not find HBHE rechits "<<inputTagHcalRecHitsHBHE_;
     LogError("PFRecHitProducerHCAL")<<err.str()<<endl;
     
     throw cms::Exception( "MissingProduct", err.str());
       }
       else {
     assert( hbheHandle.isValid() );
       }
       
       // create rechits
       typedef CaloTowerCollection::const_iterator ICT;
     
       for(ICT ict=caloTowers->begin(); ict!=caloTowers->end();ict++) {
       
     const CaloTower& ct = (*ict);
       
     //C 
     if(!caloTowerGeometry) 
       caloTowerGeometry = geoHandle->getSubdetectorGeometry(ct.id());
 
       
     // get the hadronic energy.
     
     // Mike: Just ask for the Hadronic part only now!
     // Patrick : ARGH ! While this is ok for the HCAL, this is 
     // just wrong for the HF (in which em/had are artificially 
     // separated. 
     double energy = ct.hadEnergy();
     //Auguste: Photons in HF have no hadEnergy in fastsim: -> all RecHit collections are empty with photons.
     double energyEM = ct.emEnergy(); // For HF !
     //so test the total energy to deal with the photons in  HF:
     if( (energy+energyEM) < 1e-9 ) continue;
       
     assert( ct.constituentsSize() );      
     //Mike: The DetId will be taken by the first Hadronic constituent
     //         of the tower. That is only what we need
 
     
     //get the constituents of the tower
     const std::vector<DetId>& hits = ct.constituents();
     const std::vector<DetId>& allConstituents = theTowerConstituentsMap->constituentsOf(ct.id());
 
     /*
     for(unsigned int i=0;i< hits.size();++i) {
       if(hits[i].det()==DetId::Hcal) {
         HcalDetId did = hits[i];
         if ( did.subdet()==HcalEndcap || did.subdet()==HcalForward ) { 
           //double en = hits[i].energy();
           int ieta = did.ieta();
           const CaloCellGeometry *thisCell = hcalEndcapGeometry->getGeometry(did);
           const GlobalPoint& position = thisCell->getPosition();
           if ( abs(ieta) > 27 && abs(ieta) < 33 && energy > 10. ) { 
         std::cout << "HE/HF hit " << i << " at eta = " << ieta 
               << " with CT energy = " << energy 
               << " at eta, z (hit) = " << position.eta() << " " << position.z()
               << " at eta, z (cte) = " << ct.emPosition().eta() << " " << ct.emPosition().z()
               << " at eta, z (cth) = " << ct.hadPosition().eta() << " " << ct.hadPosition().z()
               << " at eta, z (cto) = " << ct.eta() << " " << ct.vz() 
               << std::endl;
           }
         }
       }
     }
     */
     
     //Reserve the DetId we are looking for:
 
     HcalDetId detid;
     // EcalDetId edetid;
     bool foundHCALConstituent = false;
 
  
     //Loop on the calotower constituents and search for HCAL
     double dead = 0.;
     double alive = 0.;
     for(unsigned int i=0;i< hits.size();++i) {
       if(hits[i].det()==DetId::Hcal) { 
         foundHCALConstituent = true;
         detid = hits[i];
         // An HCAL tower was found: Look for dead ECAL channels in the same CaloTower.
         if ( ECAL_Compensate_ && energy > ECAL_Threshold_ ) {
           for(unsigned int j=0;j<allConstituents.size();++j) { 
         if ( allConstituents[j].det()==DetId::Ecal ) { 
           alive += 1.;
           EcalChannelStatus::const_iterator chIt = theEcalChStatus->find(allConstituents[j]);
           unsigned int dbStatus = chIt != theEcalChStatus->end() ? chIt->getStatusCode() : 0;
           if ( dbStatus > ECAL_Dead_Code_ ) dead += 1.;
         }
           }
         } 
         // Protection: tower 29 in HF is merged with tower 30. 
         // Just take the position of tower 30 in that case. 
         if ( detid.subdet() == HcalForward && abs(detid.ieta()) == 29 ) continue; 
         break;
       }
     }
 
     // In case of dead ECAL channel, rescale the HCAL energy...
     double rescaleFactor = alive > 0. ? 1. + ECAL_Compensation_*dead/alive : 1.;
       
     reco::PFRecHit* pfrh = 0;
     reco::PFRecHit* pfrhCleaned = 0;
     //---ab: need 2 rechits for the HF:
     reco::PFRecHit* pfrhHFEM = 0;
     reco::PFRecHit* pfrhHFHAD = 0;
     reco::PFRecHit* pfrhHFEMCleaned = 0;
     reco::PFRecHit* pfrhHFHADCleaned = 0;
     reco::PFRecHit* pfrhHFEMCleaned29 = 0;
     reco::PFRecHit* pfrhHFHADCleaned29 = 0;
 
     if(foundHCALConstituent)
       {
         // std::cout << ", new Energy = " << energy << std::endl;
         switch( detid.subdet() ) {
         case HcalBarrel: 
           {
         if(energy < thresh_Barrel_ ) continue;
 
         /*
         // Check the timing
         if ( energy > 5. ) { 
           for(unsigned int i=0;i< hits.size();++i) {
             if( hits[i].det() != DetId::Hcal ) continue; 
             HcalDetId theDetId = hits[i]; 
             typedef HBHERecHitCollection::const_iterator iHBHE;
             iHBHE theHit = hbheHandle->find(theDetId); 
             if ( theHit != hbheHandle->end() ) 
               std::cout << "HCAL hit : " 
                 << theDetId.ieta() << " " << theDetId.iphi() << " " 
                 << theHit->energy() << " " << theHit->time() << std::endl;
           }
         }
         */
 
 
         // if ( HCAL_Calib_ ) energy   *= std::min(max_Calib_,myPFCorr->getValues(detid)->getValue());
         //if ( rescaleFactor > 1. ) 
         // std::cout << "Barrel HCAL energy rescaled from = " << energy << " to " << energy*rescaleFactor << std::endl;
         if ( rescaleFactor > 1. ) { 
           pfrhCleaned = createHcalRecHit( detid, 
                           energy, 
                           PFLayer::HCAL_BARREL1, 
                           hcalBarrelGeometry,
                           ct.id().rawId() );
           pfrhCleaned->setRescale(rescaleFactor);
           energy *= rescaleFactor;
         }
         pfrh = createHcalRecHit( detid, 
                      energy, 
                  PFLayer::HCAL_BARREL1, 
                      hcalBarrelGeometry,
                      ct.id().rawId() );
         pfrh->setRescale(rescaleFactor);
           }
           break;
         case HcalEndcap:
           {
         if(energy < thresh_Endcap_ ) continue;
 
         /*
         // Check the timing
         if ( energy > 5. ) { 
           for(unsigned int i=0;i< hits.size();++i) {
             if( hits[i].det() != DetId::Hcal ) continue; 
             HcalDetId theDetId = hits[i]; 
             typedef HBHERecHitCollection::const_iterator iHBHE;
             iHBHE theHit = hbheHandle->find(theDetId); 
             if ( theHit != hbheHandle->end() ) 
               std::cout << "HCAL hit : " 
                 << theDetId.ieta() << " " << theDetId.iphi() << " " 
                 << theHit->energy() << " " << theHit->time() << std::endl;
           }
         }
         */
 
         // Apply tower 29 calibration
         if ( HCAL_Calib_ && abs(detid.ieta()) == 29 ) energy *= HCAL_Calib_29;
         //if ( rescaleFactor > 1. ) 
         // std::cout << "End-cap HCAL energy rescaled from = " << energy << " to " << energy*rescaleFactor << std::endl;
         if ( rescaleFactor > 1. ) { 
           pfrhCleaned = createHcalRecHit( detid, 
                           energy, 
                           PFLayer::HCAL_ENDCAP, 
                           hcalEndcapGeometry,
                           ct.id().rawId() );
           pfrhCleaned->setRescale(rescaleFactor);
           energy *= rescaleFactor;
         }
         pfrh = createHcalRecHit( detid, 
                      energy, 
                      PFLayer::HCAL_ENDCAP, 
                      hcalEndcapGeometry,
                      ct.id().rawId() );
         pfrh->setRescale(rescaleFactor);
           }
           break;
         case HcalOuter:
           {
           }
           break;
         case HcalForward:
           {
         //---ab: 2 rechits for HF:
         //double energyemHF = weight_HFem_*ct.emEnergy();
         //double energyhadHF = weight_HFhad_*ct.hadEnergy();
         double energyemHF = weight_HFem_ * energyEM;
         double energyhadHF = weight_HFhad_ * energy;
         // Some energy in the tower !
         if((energyemHF+energyhadHF) < thresh_HF_ ) continue;
 
         // Some cleaning in the HF 
         double longFibre = energyemHF + energyhadHF/2.;
         double shortFibre = energyhadHF/2.;
         int ieta = detid.ieta();
         int iphi = detid.iphi();
         HcalDetId theLongDetId (HcalForward, ieta, iphi, 1);
         HcalDetId theShortDetId (HcalForward, ieta, iphi, 2);
         typedef HFRecHitCollection::const_iterator iHF;
         iHF theLongHit = hfHandle->find(theLongDetId); 
         iHF theShortHit = hfHandle->find(theShortDetId); 
         // 
         double theLongHitEnergy = 0.;
         double theShortHitEnergy = 0.;
         bool flagShortDPG =  false; 
         bool flagLongDPG = false; 
         bool flagShortTimeDPG = false; 
         bool flagLongTimeDPG = false;
         bool flagShortPulseDPG = false;
         bool flagLongPulseDPG = false;
         //
         if ( theLongHit != hfHandle->end() ) { 
           int theLongFlag = theLongHit->flags();
           theLongHitEnergy = theLongHit->energy();
           flagLongDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId, theLongFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
           flagLongTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId, theLongFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
           flagLongPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId, theLongFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
 
           //flagLongDPG =  applyLongShortDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFLongShort)==1;
           //flagLongTimeDPG = applyTimeDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
           //flagLongPulseDPG = applyPulseDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
         }
         //
         if ( theShortHit != hfHandle->end() ) { 
           int theShortFlag = theShortHit->flags();
           theShortHitEnergy = theShortHit->energy();
           flagShortDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId, theShortFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
           flagShortTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId, theShortFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
           flagShortPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId, theShortFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
           //flagShortDPG =  applyLongShortDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFLongShort)==1;
           //flagShortTimeDPG = applyTimeDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
           //flagShortPulseDPG = applyPulseDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
         }
 
         // Then check the timing in short and long fibres in all other towers.
         if ( theShortHitEnergy > longShortFibre_Cut && 
              ( theShortHit->time() < minShortTiming_Cut ||
                theShortHit->time() > maxShortTiming_Cut || 
                flagShortTimeDPG || flagShortPulseDPG ) ) { 
           // rescaleFactor = 0. ;
           pfrhHFHADCleaned = createHcalRecHit( detid, 
                                theShortHitEnergy, 
                                PFLayer::HF_HAD, 
                                hcalEndcapGeometry,
                                ct.id().rawId() );
           pfrhHFHADCleaned->setRescale(theShortHit->time());
           /*
           std::cout << "ieta/iphi = " << ieta << " " << iphi 
                 << ", Energy em/had/long/short = " 
                 << energyemHF << " " << energyhadHF << " "
                 << theLongHitEnergy << " " << theShortHitEnergy << " " 
                 << ". Time = " << theShortHit->time() << " " 
                 << ". The short and long flags : " 
                 << flagShortDPG << " " << flagLongDPG << " "  
                 << flagShortTimeDPG << " " << flagLongTimeDPG << " "  
                 << flagShortPulseDPG << " " << flagLongPulseDPG << " "  
                 << ". Short fibres were cleaned." << std::endl;
           */
           shortFibre -= theShortHitEnergy;
           theShortHitEnergy = 0.;
         }
         
         
         if ( theLongHitEnergy > longShortFibre_Cut && 
              ( theLongHit->time() < minLongTiming_Cut ||
                theLongHit->time() > maxLongTiming_Cut  || 
                flagLongTimeDPG || flagLongPulseDPG ) ) { 
           //rescaleFactor = 0. ;
           pfrhHFEMCleaned = createHcalRecHit( detid, 
                               theLongHitEnergy, 
                               PFLayer::HF_EM, 
                               hcalEndcapGeometry,
                               ct.id().rawId() );
           pfrhHFEMCleaned->setRescale(theLongHit->time());
           /*
           std::cout << "ieta/iphi = " << ieta << " " << iphi 
                 << ", Energy em/had/long/short = " 
                 << energyemHF << " " << energyhadHF << " "
                 << theLongHitEnergy << " " << theShortHitEnergy << " " 
                 << ". Time = " << theLongHit->time() << " " 
                 << ". The short and long flags : " 
                 << flagShortDPG << " " << flagLongDPG << " "  
                 << flagShortTimeDPG << " " << flagLongTimeDPG << " "  
                 << flagShortPulseDPG << " " << flagLongPulseDPG << " "  
                 << ". Long fibres were cleaned." << std::endl;
           */
           longFibre -= theLongHitEnergy;
           theLongHitEnergy = 0.;
         }
 
         // Some energy must be in the long fibres is there is some energy in the short fibres ! 
         if ( theShortHitEnergy > shortFibre_Cut && 
              ( theLongHitEnergy/theShortHitEnergy < longFibre_Fraction || 
                flagShortDPG ) ) {
           // Check if the long-fibre hit was not cleaned already (because hot)
           // In this case don't apply the cleaning
           const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theLongDetId);
           unsigned theStatusValue = theStatus->getValue();
           int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
           // The channel is killed
           if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
             // rescaleFactor = 0. ;
             pfrhHFHADCleaned = createHcalRecHit( detid, 
                              theShortHitEnergy, 
                              PFLayer::HF_HAD, 
                              hcalEndcapGeometry,
                              ct.id().rawId() );
             pfrhHFHADCleaned->setRescale(theShortHit->time());
             /*
             std::cout << "ieta/iphi = " << ieta << " " << iphi 
                   << ", Energy em/had/long/short = " 
                   << energyemHF << " " << energyhadHF << " "
                   << theLongHitEnergy << " " << theShortHitEnergy << " " 
                   << ". Time = " << theShortHit->time() << " " 
                   << ". The status value is " << theStatusValue
                   << ". The short and long flags : " 
                   << flagShortDPG << " " << flagLongDPG << " "  
                   << flagShortTimeDPG << " " << flagLongTimeDPG << " "  
                   << flagShortPulseDPG << " " << flagLongPulseDPG << " "  
                   << ". Short fibres were cleaned." << std::endl;
             */
             shortFibre -= theShortHitEnergy;
             theShortHitEnergy = 0.;
           }
         }
 
         if ( theLongHitEnergy > longFibre_Cut && 
              ( theShortHitEnergy/theLongHitEnergy < shortFibre_Fraction || 
                flagLongDPG ) ) {
           // Check if the long-fibre hit was not cleaned already (because hot)
           // In this case don't apply the cleaning
           const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId);
           unsigned theStatusValue = theStatus->getValue();
 
           int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
           // The channel is killed
           if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
             
             //rescaleFactor = 0. ;
             pfrhHFEMCleaned = createHcalRecHit( detid, 
                               theLongHitEnergy, 
                               PFLayer::HF_EM, 
                               hcalEndcapGeometry,
                               ct.id().rawId() );
             pfrhHFEMCleaned->setRescale(theLongHit->time());
             /*
             std::cout << "ieta/iphi = " << ieta << " " << iphi 
                   << ", Energy em/had/long/short = " 
                   << energyemHF << " " << energyhadHF << " "
                   << theLongHitEnergy << " " << theShortHitEnergy << " " 
                   << ". The status value is " << theStatusValue
                   << ". Time = " << theLongHit->time() << " " 
                   << ". The short and long flags : " 
                   << flagShortDPG << " " << flagLongDPG << " "  
                   << flagShortTimeDPG << " " << flagLongTimeDPG << " "  
                   << flagShortPulseDPG << " " << flagLongPulseDPG << " "  
                   << ". Long fibres were cleaned." << std::endl;
             */
             longFibre -= theLongHitEnergy;
             theLongHitEnergy = 0.;
           }
         }
 
         // Special treatment for tower 29
         // A tower with energy only at ieta = +/- 29 is not physical -> Clean
         if ( abs(ieta) == 29 ) { 
           // rescaleFactor = 0. ;
           // Clean long fibres
           if ( theLongHitEnergy > shortFibre_Cut/2. ) { 
             pfrhHFEMCleaned29 = createHcalRecHit( detid, 
                               theLongHitEnergy, 
                               PFLayer::HF_EM, 
                               hcalEndcapGeometry,
                               ct.id().rawId() );
             pfrhHFEMCleaned29->setRescale(theLongHit->time());
             /*
             std::cout << "ieta/iphi = " << ieta << " " << iphi 
                   << ", Energy em/had/long/short = " 
                   << energyemHF << " " << energyhadHF << " "
                   << theLongHitEnergy << " " << theShortHitEnergy << " " 
                   << ". Long fibres were cleaned." << std::endl;
             */
             longFibre -= theLongHitEnergy;
             theLongHitEnergy = 0.;
           }
           // Clean short fibres
           if ( theShortHitEnergy > shortFibre_Cut/2. ) { 
             pfrhHFHADCleaned29 = createHcalRecHit( detid, 
                                theShortHitEnergy, 
                                PFLayer::HF_HAD, 
                                hcalEndcapGeometry,
                                ct.id().rawId() );
             pfrhHFHADCleaned29->setRescale(theShortHit->time());
             /*
             std::cout << "ieta/iphi = " << ieta << " " << iphi 
                   << ", Energy em/had/long/short = " 
                   << energyemHF << " " << energyhadHF << " "
                   << theLongHitEnergy << " " << theShortHitEnergy << " " 
                   << ". Short fibres were cleaned." << std::endl;
             */
             shortFibre -= theShortHitEnergy;
             theShortHitEnergy = 0.;
           }
         }
         // Check the timing of the long and short fibre rechits
         
         // First, check the timing of long and short fibre in eta = 29 if tower 30.
         else if ( abs(ieta) == 30 ) { 
           int ieta29 = ieta > 0 ? 29 : -29;
           HcalDetId theLongDetId29 (HcalForward, ieta29, iphi, 1);
           HcalDetId theShortDetId29 (HcalForward, ieta29, iphi, 2);
           iHF theLongHit29 = hfHandle->find(theLongDetId29); 
           iHF theShortHit29 = hfHandle->find(theShortDetId29); 
           // 
           double theLongHitEnergy29 = 0.;
           double theShortHitEnergy29 = 0.;
           bool flagShortDPG29 =  false; 
           bool flagLongDPG29 = false; 
           bool flagShortTimeDPG29 = false; 
           bool flagLongTimeDPG29 = false;
           bool flagShortPulseDPG29 = false;
           bool flagLongPulseDPG29 = false;
           //
           if ( theLongHit29 != hfHandle->end() ) { 
             int theLongFlag29 = theLongHit29->flags();
             theLongHitEnergy29 = theLongHit29->energy() ;
             flagLongDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
             flagLongTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
             flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
 
             //flagLongDPG29 = applyLongShortDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFLongShort)==1;
             //flagLongTimeDPG29 = applyTimeDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
             //flagLongPulseDPG29 = applyPulseDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
           }
           //
           if ( theShortHit29 != hfHandle->end() ) { 
             int theShortFlag29 = theShortHit29->flags();
             theShortHitEnergy29 = theShortHit29->energy();        
             flagShortDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
             flagShortTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
             flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
 
             //flagShortDPG29 = applyLongShortDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFLongShort)==1;
             //flagShortTimeDPG29 = applyTimeDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1;
             //flagShortPulseDPG29 = applyPulseDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFDigiTime)==1;
           }
 
           if ( theLongHitEnergy29 > longShortFibre_Cut && 
                ( theLongHit29->time() < minLongTiming_Cut ||
                  theLongHit29->time() > maxLongTiming_Cut ||
              flagLongTimeDPG29 || flagLongPulseDPG29 ) ) { 
             //rescaleFactor = 0. ;
             pfrhHFEMCleaned29 = createHcalRecHit( detid, 
                               theLongHitEnergy29, 
                               PFLayer::HF_EM, 
                               hcalEndcapGeometry,
                               ct.id().rawId() );
             pfrhHFEMCleaned29->setRescale(theLongHit29->time());
             /*
             std::cout << "ieta/iphi = " << ieta29 << " " << iphi 
                   << ", Energy em/had/long/short = " 
                   << energyemHF << " " << energyhadHF << " "
                   << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " 
                   << ". Time = " << theLongHit29->time() << " " 
                   << ". The short and long flags : " 
                   << flagShortDPG29 << " " << flagLongDPG29 << " "  
                   << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "  
                   << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "  
                   << ". Long fibres were cleaned." << std::endl;
             */
             longFibre -= theLongHitEnergy29;
             theLongHitEnergy29 = 0;
           }
 
           if ( theShortHitEnergy29 > longShortFibre_Cut && 
                ( theShortHit29->time() < minShortTiming_Cut ||
                  theShortHit29->time() > maxShortTiming_Cut ||
              flagShortTimeDPG29 || flagShortPulseDPG29 ) ) { 
             //rescaleFactor = 0. ;
             pfrhHFHADCleaned29 = createHcalRecHit( detid, 
                              theShortHitEnergy29, 
                              PFLayer::HF_HAD, 
                              hcalEndcapGeometry,
                              ct.id().rawId() );
             pfrhHFHADCleaned29->setRescale(theShortHit29->time());
             /*
             std::cout << "ieta/iphi = " << ieta29 << " " << iphi 
                   << ", Energy em/had/long/short = " 
                   << energyemHF << " " << energyhadHF << " "
                   << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " 
                   << ". Time = " << theShortHit29->time() << " " 
                   << ". The short and long flags : " 
                   << flagShortDPG29 << " " << flagLongDPG29 << " "  
                   << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "  
                   << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "  
                   << ". Short fibres were cleaned." << std::endl;
             */
             shortFibre -= theShortHitEnergy29;
             theShortHitEnergy29 = 0.;
           }
 
           // Some energy must be in the long fibres is there is some energy in the short fibres ! 
           if ( theShortHitEnergy29 > shortFibre_Cut && 
                ( theLongHitEnergy29/theShortHitEnergy29 < 2.*longFibre_Fraction || 
              flagShortDPG29 ) ) {
             // Check if the long-fibre hit was not cleaned already (because hot)
             // In this case don't apply the cleaning
             const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theLongDetId29);
             unsigned theStatusValue = theStatus->getValue();
 
             int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
             // The channel is killed
             if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
               //rescaleFactor = 0. ;
               pfrhHFHADCleaned29 = createHcalRecHit( detid, 
                                theShortHitEnergy29, 
                                PFLayer::HF_HAD, 
                                hcalEndcapGeometry,
                                ct.id().rawId() );
               pfrhHFHADCleaned29->setRescale(theShortHit29->time());
               /*
               std::cout << "ieta/iphi = " << ieta29 << " " << iphi 
                 << ", Energy em/had/long/short = " 
                 << energyemHF << " " << energyhadHF << " "
                 << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " 
                 << ". Time = " << theShortHit29->time() << " " 
                 << ". The status value is " << theStatusValue
                 << ". The short and long flags : " 
                 << flagShortDPG29 << " " << flagLongDPG29 << " "  
                 << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "  
                 << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "  
                 << ". Short fibres were cleaned." << std::endl;
               */
               shortFibre -= theShortHitEnergy29;
               theShortHitEnergy29 = 0.;
             }       
           }
           
           // Some energy must be in the short fibres is there is some energy in the long fibres ! 
           if ( theLongHitEnergy29 > longFibre_Cut && 
                ( theShortHitEnergy29/theLongHitEnergy29 < shortFibre_Fraction || 
              flagLongDPG29 ) ) {
             // Check if the long-fibre hit was not cleaned already (because hot)
             // In this case don't apply the cleaning
             const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId29);
             unsigned theStatusValue = theStatus->getValue();
             int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue);
             // The channel is killed
             if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
 
               //rescaleFactor = 0. ;
               pfrhHFEMCleaned29 = createHcalRecHit( detid, 
                                 theLongHitEnergy29, 
                                 PFLayer::HF_EM, 
                                 hcalEndcapGeometry,
                                 ct.id().rawId() );
               pfrhHFEMCleaned29->setRescale(theLongHit29->time());
               /* 
               std::cout << "ieta/iphi = " << ieta29 << " " << iphi 
                 << ", Energy em/had/long/short = " 
                 << energyemHF << " " << energyhadHF << " "
                 << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " 
                 << ". The status value is " << theStatusValue
                 << ". Time = " << theLongHit29->time() << " " 
                 << ". The short and long flags : " 
                 << flagShortDPG29 << " " << flagLongDPG29 << " "  
                 << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " "  
                 << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " "  
                 << ". Long fibres were cleaned." << std::endl;
               */
               longFibre -= theLongHitEnergy29;
               theLongHitEnergy29 = 0.;
             }
           }
 
           // Check that the energy in tower 29 is smaller than in tower 30
           // First in long fibres
           if ( theLongHitEnergy29 > std::max(theLongHitEnergy,shortFibre_Cut/2) ) { 
             pfrhHFEMCleaned29 = createHcalRecHit( detid, 
                               theLongHitEnergy29, 
                               PFLayer::HF_EM, 
                               hcalEndcapGeometry,
                               ct.id().rawId() );
             pfrhHFEMCleaned29->setRescale(theLongHit29->time());
             /*
             std::cout << "ieta/iphi = " << ieta29 << " " << iphi 
                   << ", Energy L29/S29/L30/S30 = " 
                   << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
                   << theLongHitEnergy << " " << theShortHitEnergy << " " 
                   << ". Long fibres were cleaned." << std::endl;
             */
             longFibre -= theLongHitEnergy29;
             theLongHitEnergy29 = 0.;
           }
           // Second in short fibres
           if ( theShortHitEnergy29 > std::max(theShortHitEnergy,shortFibre_Cut/2.) ) { 
             pfrhHFHADCleaned29 = createHcalRecHit( detid, 
                                theShortHitEnergy29, 
                                PFLayer::HF_HAD, 
                                hcalEndcapGeometry,
                                ct.id().rawId() );
             pfrhHFHADCleaned29->setRescale(theShortHit29->time());
             /*
             std::cout << "ieta/iphi = " << ieta << " " << iphi 
                   << ", Energy L29/S29/L30/S30 = " 
                   << theLongHitEnergy29 << " " << theShortHitEnergy29 << " "
                   << theLongHitEnergy << " " << theShortHitEnergy << " " 
                   << ". Short fibres were cleaned." << std::endl;
             */
             shortFibre -= theShortHitEnergy29;
             theShortHitEnergy29 = 0.;
           }
         }
 
 
         // Determine EM and HAD after cleaning of short and long fibres
         energyhadHF = 2.*shortFibre;
         energyemHF = longFibre - shortFibre;
 
         // The EM energy might be negative, as it amounts to Long - Short
         // In that case, put the EM "energy" in the HAD energy
         // Just to avoid systematic positive bias due to "Short" high fluctuations
         if ( energyemHF < thresh_HF_ ) { 
           energyhadHF += energyemHF;
           energyemHF = 0.;
         }
 
         // Apply HCAL calibration factors flor towers close to 29, if requested
         if ( HF_Calib_ && abs(detid.ieta()) <= 32 ) { 
           energyhadHF *= HF_Calib_29;
           energyemHF *= HF_Calib_29;
         }
                 
         // Create an EM and a HAD rechit if above threshold.
         if ( energyemHF > thresh_HF_ || energyhadHF > thresh_HF_ ) { 
           pfrhHFEM = createHcalRecHit( detid, 
                            energyemHF, 
                            PFLayer::HF_EM, 
                            hcalEndcapGeometry,
                            ct.id().rawId() );
           pfrhHFHAD = createHcalRecHit( detid, 
                         energyhadHF, 
                         PFLayer::HF_HAD, 
                         hcalEndcapGeometry,
                         ct.id().rawId() );
           pfrhHFEM->setEnergyUp(energyhadHF);
           pfrhHFHAD->setEnergyUp(energyemHF);
         }
         
           }
           break;
         default:
           LogError("PFRecHitProducerHCAL")
         <<"CaloTower constituent: unknown layer : "
         <<detid.subdet()<<endl;
         } 
 
         if(pfrh) { 
           rechits.push_back( *pfrh );
           delete pfrh;
           idSortedRecHits.insert( make_pair(ct.id().rawId(), 
                         rechits.size()-1 ) ); 
         }
         if(pfrhCleaned) { 
           rechitsCleaned.push_back( *pfrhCleaned );
           delete pfrhCleaned;
         }
         //---ab: 2 rechits for HF:     
         if(pfrhHFEM) { 
           HFEMRecHits->push_back( *pfrhHFEM );
           delete pfrhHFEM;
           idSortedRecHitsHFEM.insert( make_pair(ct.id().rawId(), 
                         HFEMRecHits->size()-1 ) ); 
         }
         if(pfrhHFHAD) { 
           HFHADRecHits->push_back( *pfrhHFHAD );
           delete pfrhHFHAD;
           idSortedRecHitsHFHAD.insert( make_pair(ct.id().rawId(), 
                         HFHADRecHits->size()-1 ) ); 
         }
         //---ab    
         if(pfrhHFEMCleaned) { 
           rechitsCleaned.push_back( *pfrhHFEMCleaned );
           delete pfrhHFEMCleaned;
         }
         if(pfrhHFHADCleaned) { 
           rechitsCleaned.push_back( *pfrhHFHADCleaned );
           delete pfrhHFHADCleaned;
         }
         if(pfrhHFEMCleaned29) { 
           rechitsCleaned.push_back( *pfrhHFEMCleaned29 );
           delete pfrhHFEMCleaned29;
         }
         if(pfrhHFHADCleaned29) { 
           rechitsCleaned.push_back( *pfrhHFHADCleaned29 );
           delete pfrhHFHADCleaned29;
         }
       }
       }
       // do navigation 
       for(unsigned i=0; i<rechits.size(); i++ ) {
     findRecHitNeighboursCT( rechits[i], 
                 idSortedRecHits, 
                 caloTowerTopology);
       }
       for(unsigned i=0; i<HFEMRecHits->size(); i++ ) {
     findRecHitNeighboursCT( (*HFEMRecHits)[i], 
                 idSortedRecHitsHFEM, 
                 caloTowerTopology);
       }
       for(unsigned i=0; i<HFHADRecHits->size(); i++ ) {
     findRecHitNeighboursCT( (*HFHADRecHits)[i], 
                 idSortedRecHitsHFHAD, 
                 caloTowerTopology);
       }
       iEvent.put( HFHADRecHits,"HFHAD" );   
       iEvent.put( HFEMRecHits,"HFEM" ); 
     }   
   }
   else if( !(inputTagHcalRecHitsHBHE_ == InputTag()) ) { 
     // clustering is not done on CaloTowers but on HCAL rechits.
        
   
     // get the hcal topology
     HcalTopology hcalTopology;
     
     // HCAL rechits 
     //    vector<edm::Handle<HBHERecHitCollection> > hcalHandles;  
     edm::Handle<HBHERecHitCollection>  hcalHandle;  
 
     
     bool found = iEvent.getByLabel(inputTagHcalRecHitsHBHE_, 
                    hcalHandle );
 
     if(!found) {
       ostringstream err;
       err<<"could not find rechits "<<inputTagHcalRecHitsHBHE_;
       LogError("PFRecHitProducerHCAL")<<err.str()<<endl;
     
       throw cms::Exception( "MissingProduct", err.str());
     }
     else {
       assert( hcalHandle.isValid() );
       
       const edm::Handle<HBHERecHitCollection>& handle = hcalHandle;
       for(unsigned irechit=0; irechit<handle->size(); irechit++) {
     const HBHERecHit& hit = (*handle)[irechit];
     
     double energy = hit.energy();
     
     reco::PFRecHit* pfrh = 0;
     
 
     const HcalDetId& detid = hit.detid();
     switch( detid.subdet() ) {
     case HcalBarrel:
       {
         if(energy < thresh_Barrel_ ) continue;
         pfrh = createHcalRecHit( detid, 
                      energy, 
                      PFLayer::HCAL_BARREL1, 
                      hcalBarrelGeometry );
       }
       break;
     case HcalEndcap:
       {
         if(energy < thresh_Endcap_ ) continue;
         pfrh = createHcalRecHit( detid, 
                      energy, 
                      PFLayer::HCAL_ENDCAP, 
                      hcalEndcapGeometry );    
       }
       break;
     case HcalForward:
       {
         if(energy < thresh_HF_ ) continue;
         pfrh = createHcalRecHit( detid, 
                      energy, 
                      PFLayer::HF_HAD, 
                      hcalEndcapGeometry );
       }
       break;
     default:
       LogError("PFRecHitProducerHCAL")
         <<"HCAL rechit: unknown layer : "<<detid.subdet()<<endl;
       continue;
     } 
 
     if(pfrh) { 
       rechits.push_back( *pfrh );
       delete pfrh;
       idSortedRecHits.insert( make_pair(detid.rawId(), 
                         rechits.size()-1 ) ); 
     }
       }
       
       
       // do navigation:
       for(unsigned i=0; i<rechits.size(); i++ ) {
     
     findRecHitNeighbours( rechits[i], idSortedRecHits, 
                   hcalTopology, 
                   *hcalBarrelGeometry, 
                   hcalTopology,
                   *hcalEndcapGeometry);
       } // loop for navigation
     }  // endif hcal rechits were found
   } // endif clustering on rechits in hcal
 }

void PFRecHitProducerHCAL::findRecHitNeighbours	(	reco::PFRecHit &	rh,
		const std::map< unsigned, unsigned > &	sortedHits,
		const CaloSubdetectorTopology &	barrelTopo,
		const CaloSubdetectorGeometry &	barrelGeom,
		const CaloSubdetectorTopology &	endcapTopo,
		const CaloSubdetectorGeometry &	endcapGeom
	)

private

find and set the neighbours to a given rechit this works for ecal, hcal, ps

Definition at line 1117 of file PFRecHitProducerHCAL.cc.

References reco::PFRecHit::add4Neighbour(), reco::PFRecHit::add8Neighbour(), cond::rpcobgas::detid, reco::PFRecHit::detId(), east, CaloNavigator< T >::east(), PFLayer::ECAL_BARREL, PFLayer::ECAL_ENDCAP, geometry, PFLayer::HCAL_BARREL1, PFLayer::HCAL_ENDCAP, PFLayer::HF_EM, PFLayer::HF_HAD, CaloNavigator< T >::home(), i, reco::PFRecHit::layer(), north, CaloNavigator< T >::north(), PFLayer::PS1, PFLayer::PS2, DetId::rawId(), south, CaloNavigator< T >::south(), west, and CaloNavigator< T >::west().

Referenced by createRecHits().

                                                   {
   
   //cout<<"------PFRecHitProducerHcaL:findRecHitNeighbours navigation value "<<navigation_HF_<<endl;
  if(navigation_HF_ == false){
     if( rh.layer() == PFLayer::HF_HAD )
       return;
     if( rh.layer() == PFLayer::HF_EM )
       return;
   } 
   DetId detid( rh.detId() );
 
   const CaloSubdetectorTopology* topology = 0;
   const CaloSubdetectorGeometry* geometry = 0;
   
   switch( rh.layer() ) {
   case PFLayer::ECAL_ENDCAP: 
     topology = &endcapTopology;
     geometry = &endcapGeometry;
     break;
   case PFLayer::ECAL_BARREL: 
     topology = &barrelTopology;
     geometry = &barrelGeometry;
     break;
   case PFLayer::HCAL_ENDCAP:
     topology = &endcapTopology;
     geometry = &endcapGeometry;
     break;
   case PFLayer::HCAL_BARREL1:
     topology = &barrelTopology;
     geometry = &barrelGeometry;
     break;
   case PFLayer::PS1:
   case PFLayer::PS2:
     topology = &barrelTopology;
     geometry = &barrelGeometry;
     break;
   default:
     assert(0);
   }
   
   assert( topology && geometry );
 
   CaloNavigator<DetId> navigator(detid, topology);
 
   DetId north = navigator.north();  
   
   DetId northeast(0);
   if( north != DetId(0) ) {
     northeast = navigator.east();  
   }
   navigator.home();
 
 
   DetId south = navigator.south();
 
   
 
   DetId southwest(0); 
   if( south != DetId(0) ) {
     southwest = navigator.west();
   }
   navigator.home();
 
 
   DetId east = navigator.east();
   DetId southeast;
   if( east != DetId(0) ) {
     southeast = navigator.south(); 
   }
   navigator.home();
   DetId west = navigator.west();
   DetId northwest;
   if( west != DetId(0) ) {   
     northwest = navigator.north();  
   }
   navigator.home();
     
   IDH i = sortedHits.find( north.rawId() );
   if(i != sortedHits.end() ) 
     rh.add4Neighbour( i->second );
   
   i = sortedHits.find( northeast.rawId() );
   if(i != sortedHits.end() ) 
     rh.add8Neighbour( i->second );
   
   i = sortedHits.find( south.rawId() );
   if(i != sortedHits.end() ) 
     rh.add4Neighbour( i->second );
     
   i = sortedHits.find( southwest.rawId() );
   if(i != sortedHits.end() ) 
     rh.add8Neighbour( i->second );
     
   i = sortedHits.find( east.rawId() );
   if(i != sortedHits.end() ) 
     rh.add4Neighbour( i->second );
     
   i = sortedHits.find( southeast.rawId() );
   if(i != sortedHits.end() ) 
     rh.add8Neighbour( i->second );
     
   i = sortedHits.find( west.rawId() );
   if(i != sortedHits.end() ) 
      rh.add4Neighbour( i->second );
    
   i = sortedHits.find( northwest.rawId() );
   if(i != sortedHits.end() ) 
     rh.add8Neighbour( i->second );
     
 
 }

void PFRecHitProducerHCAL::findRecHitNeighboursCT	(	reco::PFRecHit &	rh,
		const std::map< unsigned, unsigned > &	sortedHits,
		const CaloSubdetectorTopology &	topology
	)

private

find and set the neighbours to a given rechit this works for hcal CaloTowers. Should be possible to have a single function for all detectors

Definition at line 1237 of file PFRecHitProducerHCAL.cc.

References reco::PFRecHit::add4Neighbour(), reco::PFRecHit::add8Neighbour(), reco::PFRecHit::detId(), east, CaloSubdetectorTopology::east(), PFLayer::HF_EM, PFLayer::HF_HAD, i, reco::PFRecHit::layer(), north, CaloSubdetectorTopology::north(), DetId::rawId(), south, CaloSubdetectorTopology::south(), west, and CaloSubdetectorTopology::west().

Referenced by createRecHits().

                                             {
   //cout<<"------PFRecHitProducerHcaL:findRecHitNeighboursCT navigation value "<<navigation_HF_<<endl;
   //  cout<<"----------- rechit print out"<<endl;
   // if(( rh.layer() == PFLayer::HF_HAD )||(rh.layer() == PFLayer::HF_EM)) {  
     
   //    cout<<rh<<endl;
     //  }
   if(navigation_HF_ == false){
     if( rh.layer() == PFLayer::HF_HAD )
       return;
     if( rh.layer() == PFLayer::HF_EM )
       return;
   }
   CaloTowerDetId ctDetId( rh.detId() );
     
 
   vector<DetId> northids = topology.north(ctDetId);
   vector<DetId> westids = topology.west(ctDetId);
   vector<DetId> southids = topology.south(ctDetId);
   vector<DetId> eastids = topology.east(ctDetId);
 
 
   CaloTowerDetId badId;
 
   // all the following detids will be CaloTowerDetId
   CaloTowerDetId north;
   CaloTowerDetId northwest;
   CaloTowerDetId northwest2;
   CaloTowerDetId west;
   CaloTowerDetId west2;
   CaloTowerDetId southwest;
   CaloTowerDetId southwest2;
   CaloTowerDetId south;
   CaloTowerDetId southeast;
   CaloTowerDetId southeast2;
   CaloTowerDetId east;
   CaloTowerDetId east2;
   CaloTowerDetId northeast;
   CaloTowerDetId northeast2;
   
   // for north and south, there is no ambiguity : 1 or 0 neighbours
   
   switch( northids.size() ) {
   case 0: 
     break;
   case 1: 
     north = northids[0];
     break;
   default:
   stringstream err("PFRecHitProducerHCAL::findRecHitNeighboursCT : incorrect number of neighbours north: "); 
     err<<northids.size();
     throw( err.str() ); 
   }
 
   switch( southids.size() ) {
   case 0: 
     break;
   case 1: 
     south = southids[0];
     break;
   default:
   stringstream err("PFRecHitProducerHCAL::findRecHitNeighboursCT : incorrect number of neighbours south: "); 
     err<<southids.size();
     throw( err.str() ); 
   }
   
   // for east and west, one must take care 
   // of the pitch change in HCAL endcap.
 
   switch( eastids.size() ) {
   case 0: 
     break;
   case 1: 
     east = eastids[0];
     northeast = getNorth(east, topology);
     southeast = getSouth(east, topology);
     break;
   case 2:  
     // in this case, 0 is more on the north than 1
     east = eastids[0];
     east2 = eastids[1];
     northeast = getNorth(east, topology );
     southeast = getSouth(east2, topology);    
     northeast2 = getNorth(northeast, topology );
     southeast2 = getSouth(southeast, topology);    
     break;
   default:
   stringstream err("PFRecHitProducerHCAL::findRecHitNeighboursCT : incorrect number of neighbours eastids: "); 
     err<<eastids.size();
     throw( err.str() ); 
   }
   
   
   switch( westids.size() ) {
   case 0: 
     break;
   case 1: 
     west = westids[0];
     northwest = getNorth(west, topology);
     southwest = getSouth(west, topology);
     break;
   case 2:  
     // in this case, 0 is more on the north than 1
     west = westids[0];
     west2 = westids[1];
     northwest = getNorth(west, topology );
     southwest = getSouth(west2, topology );    
     northwest2 = getNorth(northwest, topology );
     southwest2 = getSouth(southwest, topology );    
     break;
   default:
   stringstream err("PFRecHitProducerHCAL::findRecHitNeighboursCT : incorrect number of neighbours westids: "); 
     err<< westids.size();
     throw( err.str() ); 
   }
 
 
 
 
   // find and set neighbours
     
   IDH i = sortedHits.find( north.rawId() );
   if(i != sortedHits.end() ) 
     rh.add4Neighbour( i->second );
   
   i = sortedHits.find( northeast.rawId() );
   if(i != sortedHits.end() ) 
     rh.add8Neighbour( i->second );
   
   i = sortedHits.find( northeast2.rawId() );
   if(i != sortedHits.end() ) 
     rh.add8Neighbour( i->second );
   
   i = sortedHits.find( south.rawId() );
   if(i != sortedHits.end() ) 
     rh.add4Neighbour( i->second );
     
   i = sortedHits.find( southwest.rawId() );
   if(i != sortedHits.end() ) 
     rh.add8Neighbour( i->second );
     
   i = sortedHits.find( southwest2.rawId() );
   if(i != sortedHits.end() ) 
     rh.add8Neighbour( i->second );
     
   i = sortedHits.find( east.rawId() );
   if(i != sortedHits.end() ) 
     rh.add4Neighbour( i->second );
     
   i = sortedHits.find( east2.rawId() );
   if(i != sortedHits.end() ) 
     rh.add4Neighbour( i->second );
     
   i = sortedHits.find( southeast.rawId() );
   if(i != sortedHits.end() ) 
     rh.add8Neighbour( i->second );
     
   i = sortedHits.find( southeast2.rawId() );
   if(i != sortedHits.end() ) 
     rh.add8Neighbour( i->second );
     
   i = sortedHits.find( west.rawId() );
   if(i != sortedHits.end() ) 
      rh.add4Neighbour( i->second );
    
   i = sortedHits.find( west2.rawId() );
   if(i != sortedHits.end() ) 
      rh.add4Neighbour( i->second );
    
   i = sortedHits.find( northwest.rawId() );
   if(i != sortedHits.end() ) 
     rh.add8Neighbour( i->second );
 
   i = sortedHits.find( northwest2.rawId() );
   if(i != sortedHits.end() ) 
     rh.add8Neighbour( i->second );
 
   //  cout<<"----------- rechit print out"<<endl;
   // if(( rh.layer() == PFLayer::HF_HAD )||(rh.layer() == PFLayer::HF_EM)) {  
     
   //   cout<<rh<<endl;
     //  }
 }

DetId PFRecHitProducerHCAL::getNorth	(	const DetId &	id,
		const CaloSubdetectorTopology &	topology
	)

private

Definition at line 1440 of file PFRecHitProducerHCAL.cc.

References north, and CaloSubdetectorTopology::north().

                                                             {
 
   DetId north;
   vector<DetId> nids = topology.north(id);
   if(nids.size() == 1)
     north = nids[0];
   
   return north;
 } 

DetId PFRecHitProducerHCAL::getSouth	(	const DetId &	id,
		const CaloSubdetectorTopology &	topology
	)

private

Definition at line 1426 of file PFRecHitProducerHCAL.cc.

References south, and CaloSubdetectorTopology::south().

                                                             {
 
   DetId south;
   vector<DetId> sids = topology.south(id);
   if(sids.size() == 1)
     south = sids[0];
   
   return south;
 } 