72 shortFibre_Cut = iConfig.
getParameter<
double>(
"ShortFibre_Cut");
73 longFibre_Fraction = iConfig.
getParameter<
double>(
"LongFibre_Fraction");
75 longFibre_Cut = iConfig.
getParameter<
double>(
"LongFibre_Cut");
76 shortFibre_Fraction = iConfig.
getParameter<
double>(
"ShortFibre_Fraction");
78 applyLongShortDPG_ = iConfig.
getParameter<
bool>(
"ApplyLongShortDPG");
80 longShortFibre_Cut = iConfig.
getParameter<
double>(
"LongShortFibre_Cut");
81 minShortTiming_Cut = iConfig.
getParameter<
double>(
"MinShortTiming_Cut");
82 maxShortTiming_Cut = iConfig.
getParameter<
double>(
"MaxShortTiming_Cut");
83 minLongTiming_Cut = iConfig.
getParameter<
double>(
"MinLongTiming_Cut");
84 maxLongTiming_Cut = iConfig.
getParameter<
double>(
"MaxLongTiming_Cut");
86 applyTimeDPG_ = iConfig.
getParameter<
bool>(
"ApplyTimeDPG");
87 applyPulseDPG_ = iConfig.
getParameter<
bool>(
"ApplyPulseDPG");
88 HcalMaxAllowedHFLongShortSev_ = iConfig.
getParameter<
int>(
"HcalMaxAllowedHFLongShortSev");
89 HcalMaxAllowedHFDigiTimeSev_ = iConfig.
getParameter<
int>(
"HcalMaxAllowedHFDigiTimeSev");
90 HcalMaxAllowedHFInTimeWindowSev_ = iConfig.
getParameter<
int>(
"HcalMaxAllowedHFInTimeWindowSev");
91 HcalMaxAllowedChannelStatusSev_ = iConfig.
getParameter<
int>(
"HcalMaxAllowedChannelStatusSev");
93 ECAL_Compensate_ = iConfig.
getParameter<
bool>(
"ECAL_Compensate");
94 ECAL_Threshold_ = iConfig.
getParameter<
double>(
"ECAL_Threshold");
95 ECAL_Compensation_ = iConfig.
getParameter<
double>(
"ECAL_Compensation");
96 ECAL_Dead_Code_ = iConfig.
getParameter<
unsigned int>(
"ECAL_Dead_Code");
115 produces<reco::PFRecHitCollection>();
116 produces<reco::PFRecHitCollection>(
"Cleaned");
117 produces<reco::PFRecHitCollection>(
"HFHAD").setBranchAlias(
"HFHADRecHits");
118 produces<reco::PFRecHitCollection>(
"HFEM").setBranchAlias(
"HFEMRecHits");
127 navigator_->beginEvent(iSetup);
148 theHcalTopology = hcalTopology.
product();
150 auto rechits = std::make_unique<std::vector<reco::PFRecHit>>();
151 auto rechitsCleaned = std::make_unique<std::vector<reco::PFRecHit>>();
152 auto HFHADRecHits = std::make_unique<std::vector<reco::PFRecHit>>();
153 auto HFEMRecHits = std::make_unique<std::vector<reco::PFRecHit>>();
165 for(
auto const & ct : *caloTowers) {
172 double energy = ct.hadEnergy();
174 double energyEM = ct.emEnergy();
176 if( (energy+energyEM) < 1
e-9 )
continue;
179 const std::vector<DetId>&
hits = ct.constituents();
180 const std::vector<DetId>& allConstituents = theTowerConstituentsMap->constituentsOf(ct.id());
182 bool foundHCALConstituent =
false;
186 for(
auto hit : hits) {
188 foundHCALConstituent =
true;
190 if (theHcalTopology->getMergePositionFlag() &&
192 detid = theHcalTopology->idFront(detid);
195 if ( ECAL_Compensate_ && energy > ECAL_Threshold_ ) {
196 for(
auto allConstituent : allConstituents) {
199 auto chIt = theEcalChStatus->find(allConstituent);
200 unsigned int dbStatus = chIt != theEcalChStatus->end() ? chIt->getStatusCode() : 0;
201 if ( dbStatus > ECAL_Dead_Code_ ) dead += 1.;
213 double rescaleFactor = alive > 0. ? 1. + ECAL_Compensation_*dead/alive : 1.;
225 if(foundHCALConstituent)
228 switch( detid.
subdet() ) {
231 if(energy < thresh_Barrel_ )
continue;
232 if ( rescaleFactor > 1. ) {
233 pfrhCleaned = createHcalRecHit( detid,
238 pfrhCleaned->
setTime(rescaleFactor);
239 energy *= rescaleFactor;
241 pfrh = createHcalRecHit( detid,
251 if(energy < thresh_Endcap_ )
continue;
254 if ( rescaleFactor > 1. ) {
255 pfrhCleaned = createHcalRecHit( detid,
260 pfrhCleaned->
setTime(rescaleFactor);
261 energy *= rescaleFactor;
263 pfrh = createHcalRecHit( detid,
280 double energyemHF = weight_HFem_ * energyEM;
281 double energyhadHF = weight_HFhad_ * energy;
283 if((energyemHF+energyhadHF) < thresh_HF_ )
continue;
286 double longFibre = energyemHF + energyhadHF/2.;
287 double shortFibre = energyhadHF/2.;
288 int ieta = detid.
ieta();
289 int iphi = detid.
iphi();
293 auto theLongHit = hfHandle->
find(theLongDetId);
294 auto theShortHit = hfHandle->
find(theShortDetId);
296 double theLongHitEnergy = 0.;
297 double theShortHitEnergy = 0.;
298 bool flagShortDPG =
false;
299 bool flagLongDPG =
false;
300 bool flagShortTimeDPG =
false;
301 bool flagLongTimeDPG =
false;
302 bool flagShortPulseDPG =
false;
303 bool flagLongPulseDPG =
false;
305 if ( theLongHit != hfHandle->
end() ) {
306 int theLongFlag = theLongHit->flags();
307 theLongHitEnergy = theLongHit->energy();
308 flagLongDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId, theLongFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
309 flagLongTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId, theLongFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
310 flagLongPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId, theLongFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
317 if ( theShortHit != hfHandle->
end() ) {
318 int theShortFlag = theShortHit->flags();
319 theShortHitEnergy = theShortHit->energy();
320 flagShortDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId, theShortFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
321 flagShortTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId, theShortFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
322 flagShortPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId, theShortFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
329 if ( theShortHitEnergy > longShortFibre_Cut &&
330 ( theShortHit->time() < minShortTiming_Cut ||
331 theShortHit->time() > maxShortTiming_Cut ||
332 flagShortTimeDPG || flagShortPulseDPG ) ) {
334 pfrhHFHADCleaned = createHcalRecHit( detid,
339 pfrhHFHADCleaned->
setTime(theShortHit->time());
352 shortFibre -= theShortHitEnergy;
353 theShortHitEnergy = 0.;
357 if ( theLongHitEnergy > longShortFibre_Cut &&
358 ( theLongHit->time() < minLongTiming_Cut ||
359 theLongHit->time() > maxLongTiming_Cut ||
360 flagLongTimeDPG || flagLongPulseDPG ) ) {
362 pfrhHFEMCleaned = createHcalRecHit( detid,
367 pfrhHFEMCleaned->
setTime(theLongHit->time());
380 longFibre -= theLongHitEnergy;
381 theLongHitEnergy = 0.;
385 if ( theShortHitEnergy > shortFibre_Cut &&
386 ( theLongHitEnergy/theShortHitEnergy < longFibre_Fraction ||
391 unsigned theStatusValue = theStatus->
getValue();
392 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
395 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
397 pfrhHFHADCleaned = createHcalRecHit( detid,
402 pfrhHFHADCleaned->
setTime(theShortHit->time());
416 shortFibre -= theShortHitEnergy;
417 theShortHitEnergy = 0.;
421 if ( theLongHitEnergy > longFibre_Cut &&
422 ( theShortHitEnergy/theLongHitEnergy < shortFibre_Fraction ||
427 unsigned theStatusValue = theStatus->
getValue();
429 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
432 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
435 pfrhHFEMCleaned = createHcalRecHit( detid,
440 pfrhHFEMCleaned->
setTime(theLongHit->time());
454 longFibre -= theLongHitEnergy;
455 theLongHitEnergy = 0.;
461 if (
abs(ieta) == 29 ) {
464 if ( theLongHitEnergy > shortFibre_Cut/2. ) {
465 pfrhHFEMCleaned29 = createHcalRecHit( detid,
470 pfrhHFEMCleaned29->
setTime(theLongHit->time());
478 longFibre -= theLongHitEnergy;
479 theLongHitEnergy = 0.;
482 if ( theShortHitEnergy > shortFibre_Cut/2. ) {
483 pfrhHFHADCleaned29 = createHcalRecHit( detid,
488 pfrhHFHADCleaned29->
setTime(theShortHit->time());
496 shortFibre -= theShortHitEnergy;
497 theShortHitEnergy = 0.;
503 else if (
abs(ieta) == 30 ) {
504 int ieta29 = ieta > 0 ? 29 : -29;
507 auto theLongHit29 = hfHandle->
find(theLongDetId29);
508 auto theShortHit29 = hfHandle->
find(theShortDetId29);
510 double theLongHitEnergy29 = 0.;
511 double theShortHitEnergy29 = 0.;
512 bool flagShortDPG29 =
false;
513 bool flagLongDPG29 =
false;
514 bool flagShortTimeDPG29 =
false;
515 bool flagLongTimeDPG29 =
false;
516 bool flagShortPulseDPG29 =
false;
517 bool flagLongPulseDPG29 =
false;
519 if ( theLongHit29 != hfHandle->
end() ) {
520 int theLongFlag29 = theLongHit29->flags();
521 theLongHitEnergy29 = theLongHit29->energy() ;
522 flagLongDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
523 flagLongTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
524 flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
531 if ( theShortHit29 != hfHandle->
end() ) {
532 int theShortFlag29 = theShortHit29->flags();
533 theShortHitEnergy29 = theShortHit29->energy();
534 flagShortDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
535 flagShortTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
536 flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
543 if ( theLongHitEnergy29 > longShortFibre_Cut &&
544 ( theLongHit29->time() < minLongTiming_Cut ||
545 theLongHit29->time() > maxLongTiming_Cut ||
546 flagLongTimeDPG29 || flagLongPulseDPG29 ) ) {
548 pfrhHFEMCleaned29 = createHcalRecHit( detid,
553 pfrhHFEMCleaned29->
setTime(theLongHit29->time());
566 longFibre -= theLongHitEnergy29;
567 theLongHitEnergy29 = 0;
570 if ( theShortHitEnergy29 > longShortFibre_Cut &&
571 ( theShortHit29->time() < minShortTiming_Cut ||
572 theShortHit29->time() > maxShortTiming_Cut ||
573 flagShortTimeDPG29 || flagShortPulseDPG29 ) ) {
575 pfrhHFHADCleaned29 = createHcalRecHit( detid,
580 pfrhHFHADCleaned29->
setTime(theShortHit29->time());
593 shortFibre -= theShortHitEnergy29;
594 theShortHitEnergy29 = 0;
598 if ( theShortHitEnergy29 > shortFibre_Cut &&
599 ( theLongHitEnergy29/theShortHitEnergy29 < 2.*longFibre_Fraction ||
604 unsigned theStatusValue = theStatus->
getValue();
606 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
609 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
611 pfrhHFHADCleaned29 = createHcalRecHit( detid,
616 pfrhHFHADCleaned29->
setTime(theShortHit29->time());
630 shortFibre -= theShortHitEnergy29;
631 theShortHitEnergy29 = 0.;
636 if ( theLongHitEnergy29 > longFibre_Cut &&
637 ( theShortHitEnergy29/theLongHitEnergy29 < shortFibre_Fraction ||
641 const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId29);
642 unsigned theStatusValue = theStatus->
getValue();
643 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
646 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
649 pfrhHFEMCleaned29 = createHcalRecHit( detid,
654 pfrhHFEMCleaned29->
setTime(theLongHit29->time());
668 longFibre -= theLongHitEnergy29;
669 theLongHitEnergy29 = 0.;
676 if ( theLongHitEnergy29 >
std::max(theLongHitEnergy,shortFibre_Cut/2) ) {
677 pfrhHFEMCleaned29 = createHcalRecHit( detid,
682 pfrhHFEMCleaned29->
setTime(theLongHit29->time());
690 longFibre -= theLongHitEnergy29;
691 theLongHitEnergy29 = 0.;
694 if ( theShortHitEnergy29 >
std::max(theShortHitEnergy,shortFibre_Cut/2.) ) {
695 pfrhHFHADCleaned29 = createHcalRecHit( detid,
700 pfrhHFHADCleaned29->
setTime(theShortHit29->time());
708 shortFibre -= theShortHitEnergy29;
709 theShortHitEnergy29 = 0.;
715 energyhadHF = 2.*shortFibre;
716 energyemHF = longFibre - shortFibre;
721 if ( energyemHF < thresh_HF_ ) {
722 energyhadHF += energyemHF;
727 if ( HF_Calib_ &&
abs(detid.
ieta()) <= 32 ) {
733 if ( energyemHF > thresh_HF_ || energyhadHF > thresh_HF_ ) {
734 pfrhHFEM = createHcalRecHit( detid,
739 pfrhHFHAD = createHcalRecHit( detid,
751 <<
"CaloTower constituent: unknown layer : " 761 rechitsCleaned->push_back( *pfrhCleaned );
765 HFEMRecHits->push_back( *pfrhHFEM );
769 HFHADRecHits->push_back( *pfrhHFHAD );
772 if(pfrhHFEMCleaned) {
773 rechitsCleaned->push_back( *pfrhHFEMCleaned );
774 delete pfrhHFEMCleaned;
776 if(pfrhHFHADCleaned) {
777 rechitsCleaned->push_back( *pfrhHFHADCleaned );
778 delete pfrhHFHADCleaned;
780 if(pfrhHFEMCleaned29) {
781 rechitsCleaned->push_back( *pfrhHFEMCleaned29 );
782 delete pfrhHFEMCleaned29;
784 if(pfrhHFHADCleaned29) {
785 rechitsCleaned->push_back( *pfrhHFHADCleaned29 );
786 delete pfrhHFHADCleaned29;
798 for(
unsigned int i=0;
i<
rechits->size();++
i) {
802 if (navigation_HF_) {
808 for(
unsigned int i=0;
i<HFEMRecHits->size();++
i) {
809 navigator_->associateNeighbours(HFEMRecHits->at(
i),HFEMRecHits,refProdEM);
816 for(
unsigned int i=0;
i<HFHADRecHits->size();++
i) {
817 navigator_->associateNeighbours(HFHADRecHits->at(
i),HFHADRecHits,refProdHAD);
839 theHcalChStatus = hcalChStatus.
product();
844 theEcalChStatus = ecalChStatus.
product();
848 theTowerConstituentsMap = cttopo.
product();
859 std::shared_ptr<const CaloCellGeometry> thisCell = geom->
getGeometry(detid);
862 <<
"warning detid "<<detid.
rawId()<<
" not found in layer " 871 auto zp =
dynamic_cast<IdealZPrism const*
>(thisCell.get());
873 thisCell = zp->forPF();
T getParameter(std::string const &) const
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
void beginLuminosityBlock(const edm::LuminosityBlock &lumi, const edm::EventSetup &es) override
HcalSubdetector subdet() const
get the subdetector
bool getByToken(EDGetToken token, Handle< PROD > &result) const
PFCTRecHitProducer(const edm::ParameterSet &)
constexpr uint32_t rawId() const
get the raw id
std::vector< CaloTower >::const_iterator const_iterator
~PFCTRecHitProducer() override
std::vector< PFRecHit > PFRecHitCollection
collection of PFRecHit objects
reco::PFRecHit * createHcalRecHit(const DetId &detid, double energy, PFLayer::Layer layer, const CaloSubdetectorGeometry *geom, const CaloTowerDetId &newDetId)
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
int ieta() const
get the cell ieta
Abs< T >::type abs(const T &t)
void setTime(double time)
RefProd< PROD > getRefBeforePut()
const_iterator end() const
int iphi() const
get the cell iphi
virtual std::shared_ptr< const CaloCellGeometry > getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
int getSeverityLevel(const DetId &myid, const uint32_t &myflag, const uint32_t &mystatus) const
iterator find(key_type k)
uint32_t getValue() const
T const * product() const
void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) override
T get(const Candidate &c)