|
|
|
Producer for particle flow rechits (PFRecHit) More...
#include <PFRecHitProducerPS.h>
Public Member Functions | |
| PFRecHitProducerPS (const edm::ParameterSet &) | |
| ~PFRecHitProducerPS () | |
Private Member Functions | |
| 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) |
Private Attributes | |
| edm::InputTag | inputTagEcalRecHitsES_ |
Producer for particle flow rechits (PFRecHit)
Definition at line 32 of file PFRecHitProducerPS.h.
| PFRecHitProducerPS::PFRecHitProducerPS | ( | const edm::ParameterSet & | iConfig | ) | [explicit] |
Definition at line 32 of file PFRecHitProducerPS.cc.
References edm::ParameterSet::getParameter(), and inputTagEcalRecHitsES_.
: PFRecHitProducer(iConfig) { // access to the collections of rechits inputTagEcalRecHitsES_ = iConfig.getParameter<InputTag>("ecalRecHitsES"); }
| PFRecHitProducerPS::~PFRecHitProducerPS | ( | ) |
Definition at line 45 of file PFRecHitProducerPS.cc.
{}
| void PFRecHitProducerPS::createRecHits | ( | std::vector< reco::PFRecHit > & | rechits, |
| std::vector< reco::PFRecHit > & | rechitsCleaned, | ||
| edm::Event & | , | ||
| const edm::EventSetup & | |||
| ) | [private, virtual] |
gets PS rechits, translate them to PFRecHits, which are stored in the rechits vector
if ( !theStatusValue )
if ( !theStatusValue )
if ( !theStatusValue )
if ( !theStatusValue )
Implements PFRecHitProducer.
Definition at line 115 of file PFHCALDualTimeRecHitProducer.cc.
References abs, PFHCALDualTimeRecHitProducer::applyLongShortDPG_, PFHCALDualTimeRecHitProducer::applyPulseDPG_, PFHCALDualTimeRecHitProducer::applyTimeDPG_, HiRecoJets_cff::caloTowers, CaloTower::constituents(), CaloTowerConstituentsMap::constituentsOf(), CaloTower::constituentsSize(), PFHCALDualTimeRecHitProducer::createHcalRecHit(), HcalDetId::depth(), cond::rpcobgas::detid, CaloRecHit::detid(), alignCSCRings::e, DetId::Ecal, PFHCALDualTimeRecHitProducer::ECAL_Compensate_, PFHCALDualTimeRecHitProducer::ECAL_Dead_Code_, PFHCALDualTimeRecHitProducer::ECAL_Threshold_, CaloTower::emEnergy(), EcalCondObjectContainer< T >::end(), CaloRecHit::energy(), relval_parameters_module::energy, Exception, EcalCondObjectContainer< T >::find(), PFHCALDualTimeRecHitProducer::findRecHitNeighbours(), PFHCALDualTimeRecHitProducer::findRecHitNeighboursCT(), newFWLiteAna::found, edm::EventSetup::get(), edm::Event::getByLabel(), HcalChannelStatus::getValue(), HcalCondObjectContainer< Item >::getValues(), CaloTower::hadEnergy(), patZpeak::handle, DetId::Hcal, PFLayer::HCAL_BARREL1, PFHCALDualTimeRecHitProducer::HCAL_Calib_, PFHCALDualTimeRecHitProducer::HCAL_Calib_29, PFLayer::HCAL_ENDCAP, HcalBarrel, HcalEndcap, HcalForward, HcalOuter, PFHCALDualTimeRecHitProducer::HF_Calib_, PFHCALDualTimeRecHitProducer::HF_Calib_29, PFLayer::HF_EM, PFLayer::HF_HAD, HcalCaloFlagLabels::HFDigiTime, HcalCaloFlagLabels::HFInTimeWindow, HcalCaloFlagLabels::HFLongShort, i, CaloTower::id(), HcalDetId::ieta(), PFHCALDualTimeRecHitProducer::inputTagCaloTowers_, PFHCALDualTimeRecHitProducer::inputTagHcalRecHitsHBHE_, PFHCALDualTimeRecHitProducer::inputTagHcalRecHitsHF_, HcalDetId::iphi(), edm::HandleBase::isValid(), j, PFHCALDualTimeRecHitProducer::longFibre_Cut, PFHCALDualTimeRecHitProducer::longFibre_Fraction, PFHCALDualTimeRecHitProducer::longShortFibre_Cut, max(), PFHCALDualTimeRecHitProducer::maxLongTiming_Cut, PFHCALDualTimeRecHitProducer::maxShortTiming_Cut, PFHCALDualTimeRecHitProducer::minLongTiming_Cut, PFHCALDualTimeRecHitProducer::minShortTiming_Cut, edm::Event::put(), DetId::rawId(), reco::PFRecHit::setEnergyUp(), reco::PFRecHit::setRescale(), PFHCALDualTimeRecHitProducer::shortFibre_Cut, PFHCALDualTimeRecHitProducer::shortFibre_Fraction, HcalDetId::subdet(), PFRecHitProducer::theEcalChStatus, PFRecHitProducer::theHcalChStatus, PFRecHitProducer::theTowerConstituentsMap, PFRecHitProducer::thresh_Barrel_, PFRecHitProducer::thresh_Endcap_, PFHCALDualTimeRecHitProducer::thresh_HF_, CaloRecHit::time(), PFHCALDualTimeRecHitProducer::weight_HFem_, and PFHCALDualTimeRecHitProducer::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 the hcal topology
edm::ESHandle<HcalTopology> hcalTopology;
iSetup.get<IdealGeometryRecord>().get( hcalTopology );
//--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("PFHCALDualTimeRecHitProducer")<<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("PFHCALDualTimeRecHitProducer")<<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("PFHCALDualTimeRecHitProducer")<<err.str()<<endl;
throw cms::Exception( "MissingProduct", err.str());
}
else {
assert( hbheHandle.isValid() );
}
for(unsigned irechit=0; irechit<hbheHandle->size(); irechit++) {
const HBHERecHit& hit = (*hbheHandle)[irechit];
double hitenergy = hit.energy();
double hittime = hit.time();
reco::PFRecHit* pfrh = 0;
reco::PFRecHit* pfrhCleaned = 0;
const HcalDetId& detid = hit.detid();
switch( detid.subdet() ) {
case HcalBarrel:
{
if(hitenergy < thresh_Barrel_ ) continue;
if(detid.depth()==1) {
hittime -= 48.9580/(2.16078+hitenergy);
} else if(detid.depth()==2) {
hittime -= 34.2860/(1.23746+hitenergy);
} else if(detid.depth()==3) {
hittime -= 38.6872/(1.48051+hitenergy);
}
// time window for signal=4
if( (detid.depth()==1 && hittime>-20 && hittime<5)
|| (detid.depth()==2 && hittime>-17 && hittime<8)
|| (detid.depth()==3 && hittime>-15 && hittime<10)
) {
pfrh = createHcalRecHit( detid,
hitenergy,
PFLayer::HCAL_BARREL1,
hcalBarrelGeometry );
pfrh->setRescale(hittime);
}
}
break;
case HcalEndcap:
{
if(hitenergy < thresh_Endcap_ ) continue;
// Apply tower 29 calibration
if ( HCAL_Calib_ && abs(detid.ieta()) == 29 ) hitenergy *= HCAL_Calib_29;
if(detid.depth()==1) {
hittime -= 60.8050/(3.07285+hitenergy);
} else if(detid.depth()==2) {
hittime -= 47.1677/(2.06485+hitenergy);
} else if(detid.depth()==3) {
hittime -= 37.1941/(1.53790+hitenergy);
} else if(detid.depth()==4) {
hittime -= 42.9898/(1.92969+hitenergy);
} else if(detid.depth()==5) {
hittime -= 48.3157/(2.29903+hitenergy);
}
// time window for signal=4
if( (detid.depth()==1 && hittime>-20 && hittime<5)
|| (detid.depth()==2 && hittime>-19 && hittime<6)
|| (detid.depth()==3 && hittime>-18 && hittime<7)
|| (detid.depth()==4 && hittime>-17 && hittime<8)
|| (detid.depth()==5 && hittime>-15 && hittime<10)
) {
pfrh = createHcalRecHit( detid,
hitenergy,
PFLayer::HCAL_ENDCAP,
hcalEndcapGeometry );
pfrh->setRescale(hittime);
}
}
break;
default:
LogError("PFHCALDualTimeRecHitProducer")
<<"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 ) );
}
if(pfrhCleaned) {
rechitsCleaned.push_back( *pfrhCleaned );
delete pfrhCleaned;
}
}
// 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() ) {
theLongHitEnergy = theLongHit->energy();
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() ) {
theShortHitEnergy = theShortHit->energy();
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();
// The channel is killed
if ( !theStatusValue ) {
// 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();
// The channel is killed
if ( !theStatusValue ) {
//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() ) {
theLongHitEnergy29 = theLongHit29->energy() ;
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() ) {
theShortHitEnergy29 = theShortHit29->energy();
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();
// The channel is killed
if ( !theStatusValue ) {
//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();
// The channel is killed
if ( !theStatusValue ) {
//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("PFHCALDualTimeRecHitProducer")
<<"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<rechits.size(); i++ ) {
findRecHitNeighbours( rechits[i], idSortedRecHits,
*hcalTopology,
*hcalBarrelGeometry,
*hcalTopology,
*hcalEndcapGeometry);
} // loop for navigation
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.
// 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("PFHCALDualTimeRecHitProducer")<<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("PFHCALDualTimeRecHitProducer")
<<"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 PFRecHitProducerPS::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 COLIN remonter cette fonction dans la classe de base
Definition at line 1181 of file PFHCALDualTimeRecHitProducer.cc.
References reco::PFRecHit::add4Neighbour(), reco::PFRecHit::add8Neighbour(), reco::PFRecHit::detId(), cond::rpcobgas::detid, east, PFLayer::ECAL_BARREL, PFLayer::ECAL_ENDCAP, geometry, PFLayer::HCAL_BARREL1, PFLayer::HCAL_ENDCAP, PFLayer::HF_EM, PFLayer::HF_HAD, i, reco::PFRecHit::layer(), north, PFLayer::PS1, PFLayer::PS2, DetId::rawId(), south, and west.
{
//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;
// const CaloSubdetectorGeometry* othergeometry = 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;
// othergeometry = &barrelGeometry;
break;
case PFLayer::HCAL_BARREL1:
topology = &barrelTopology;
geometry = &barrelGeometry;
// othergeometry = &endcapGeometry;
break;
case PFLayer::PS1:
case PFLayer::PS2:
topology = &barrelTopology;
geometry = &barrelGeometry;
// othergeometry = &endcapGeometry;
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 );
}
Definition at line 63 of file PFRecHitProducerPS.h.
Referenced by PFRecHitProducerPS().
1.7.3