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(ICT ict=caloTowers->begin(); ict!=caloTowers->end();ict++) {
179 if( (energy+energyEM) < 1
e-9 )
continue;
183 const std::vector<DetId>& allConstituents = theTowerConstituentsMap->constituentsOf(ct.
id());
185 bool foundHCALConstituent =
false;
189 for(
unsigned int i=0;
i< hits.size();++
i) {
191 foundHCALConstituent =
true;
193 if (theHcalTopology->withSpecialRBXHBHE() &&
195 detid = theHcalTopology->idFront(detid);
198 if ( ECAL_Compensate_ && energy > ECAL_Threshold_ ) {
199 for(
unsigned int j=0;j<allConstituents.size();++j) {
203 unsigned int dbStatus = chIt != theEcalChStatus->end() ? chIt->getStatusCode() : 0;
204 if ( dbStatus > ECAL_Dead_Code_ ) dead += 1.;
216 double rescaleFactor = alive > 0. ? 1. + ECAL_Compensation_*dead/alive : 1.;
228 if(foundHCALConstituent)
231 switch( detid.
subdet() ) {
234 if(energy < thresh_Barrel_ )
continue;
235 if ( rescaleFactor > 1. ) {
236 pfrhCleaned = createHcalRecHit( detid,
241 pfrhCleaned->
setTime(rescaleFactor);
242 energy *= rescaleFactor;
244 pfrh = createHcalRecHit( detid,
254 if(energy < thresh_Endcap_ )
continue;
257 if ( rescaleFactor > 1. ) {
258 pfrhCleaned = createHcalRecHit( detid,
263 pfrhCleaned->
setTime(rescaleFactor);
264 energy *= rescaleFactor;
266 pfrh = createHcalRecHit( detid,
283 double energyemHF = weight_HFem_ * energyEM;
284 double energyhadHF = weight_HFhad_ * energy;
286 if((energyemHF+energyhadHF) < thresh_HF_ )
continue;
289 double longFibre = energyemHF + energyhadHF/2.;
290 double shortFibre = energyhadHF/2.;
291 int ieta = detid.
ieta();
292 int iphi = detid.
iphi();
296 iHF theLongHit = hfHandle->
find(theLongDetId);
297 iHF theShortHit = hfHandle->
find(theShortDetId);
299 double theLongHitEnergy = 0.;
300 double theShortHitEnergy = 0.;
301 bool flagShortDPG =
false;
302 bool flagLongDPG =
false;
303 bool flagShortTimeDPG =
false;
304 bool flagLongTimeDPG =
false;
305 bool flagShortPulseDPG =
false;
306 bool flagLongPulseDPG =
false;
308 if ( theLongHit != hfHandle->
end() ) {
309 int theLongFlag = theLongHit->flags();
310 theLongHitEnergy = theLongHit->energy();
311 flagLongDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId, theLongFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
312 flagLongTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId, theLongFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
313 flagLongPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId, theLongFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
320 if ( theShortHit != hfHandle->
end() ) {
321 int theShortFlag = theShortHit->flags();
322 theShortHitEnergy = theShortHit->energy();
323 flagShortDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId, theShortFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
324 flagShortTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId, theShortFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
325 flagShortPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId, theShortFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
332 if ( theShortHitEnergy > longShortFibre_Cut &&
333 ( theShortHit->time() < minShortTiming_Cut ||
334 theShortHit->time() > maxShortTiming_Cut ||
335 flagShortTimeDPG || flagShortPulseDPG ) ) {
337 pfrhHFHADCleaned = createHcalRecHit( detid,
342 pfrhHFHADCleaned->
setTime(theShortHit->time());
355 shortFibre -= theShortHitEnergy;
356 theShortHitEnergy = 0.;
360 if ( theLongHitEnergy > longShortFibre_Cut &&
361 ( theLongHit->time() < minLongTiming_Cut ||
362 theLongHit->time() > maxLongTiming_Cut ||
363 flagLongTimeDPG || flagLongPulseDPG ) ) {
365 pfrhHFEMCleaned = createHcalRecHit( detid,
370 pfrhHFEMCleaned->
setTime(theLongHit->time());
383 longFibre -= theLongHitEnergy;
384 theLongHitEnergy = 0.;
388 if ( theShortHitEnergy > shortFibre_Cut &&
389 ( theLongHitEnergy/theShortHitEnergy < longFibre_Fraction ||
394 unsigned theStatusValue = theStatus->
getValue();
395 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
398 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
400 pfrhHFHADCleaned = createHcalRecHit( detid,
405 pfrhHFHADCleaned->
setTime(theShortHit->time());
419 shortFibre -= theShortHitEnergy;
420 theShortHitEnergy = 0.;
424 if ( theLongHitEnergy > longFibre_Cut &&
425 ( theShortHitEnergy/theLongHitEnergy < shortFibre_Fraction ||
430 unsigned theStatusValue = theStatus->
getValue();
432 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
435 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
438 pfrhHFEMCleaned = createHcalRecHit( detid,
443 pfrhHFEMCleaned->
setTime(theLongHit->time());
457 longFibre -= theLongHitEnergy;
458 theLongHitEnergy = 0.;
464 if (
abs(ieta) == 29 ) {
467 if ( theLongHitEnergy > shortFibre_Cut/2. ) {
468 pfrhHFEMCleaned29 = createHcalRecHit( detid,
473 pfrhHFEMCleaned29->
setTime(theLongHit->time());
481 longFibre -= theLongHitEnergy;
482 theLongHitEnergy = 0.;
485 if ( theShortHitEnergy > shortFibre_Cut/2. ) {
486 pfrhHFHADCleaned29 = createHcalRecHit( detid,
491 pfrhHFHADCleaned29->
setTime(theShortHit->time());
499 shortFibre -= theShortHitEnergy;
500 theShortHitEnergy = 0.;
506 else if (
abs(ieta) == 30 ) {
507 int ieta29 = ieta > 0 ? 29 : -29;
510 iHF theLongHit29 = hfHandle->
find(theLongDetId29);
511 iHF theShortHit29 = hfHandle->
find(theShortDetId29);
513 double theLongHitEnergy29 = 0.;
514 double theShortHitEnergy29 = 0.;
515 bool flagShortDPG29 =
false;
516 bool flagLongDPG29 =
false;
517 bool flagShortTimeDPG29 =
false;
518 bool flagLongTimeDPG29 =
false;
519 bool flagShortPulseDPG29 =
false;
520 bool flagLongPulseDPG29 =
false;
522 if ( theLongHit29 != hfHandle->
end() ) {
523 int theLongFlag29 = theLongHit29->flags();
524 theLongHitEnergy29 = theLongHit29->energy() ;
525 flagLongDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
526 flagLongTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
527 flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
534 if ( theShortHit29 != hfHandle->
end() ) {
535 int theShortFlag29 = theShortHit29->flags();
536 theShortHitEnergy29 = theShortHit29->energy();
537 flagShortDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
538 flagShortTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
539 flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
546 if ( theLongHitEnergy29 > longShortFibre_Cut &&
547 ( theLongHit29->time() < minLongTiming_Cut ||
548 theLongHit29->time() > maxLongTiming_Cut ||
549 flagLongTimeDPG29 || flagLongPulseDPG29 ) ) {
551 pfrhHFEMCleaned29 = createHcalRecHit( detid,
556 pfrhHFEMCleaned29->
setTime(theLongHit29->time());
569 longFibre -= theLongHitEnergy29;
570 theLongHitEnergy29 = 0;
573 if ( theShortHitEnergy29 > longShortFibre_Cut &&
574 ( theShortHit29->time() < minShortTiming_Cut ||
575 theShortHit29->time() > maxShortTiming_Cut ||
576 flagShortTimeDPG29 || flagShortPulseDPG29 ) ) {
578 pfrhHFHADCleaned29 = createHcalRecHit( detid,
583 pfrhHFHADCleaned29->
setTime(theShortHit29->time());
596 shortFibre -= theShortHitEnergy29;
597 theShortHitEnergy29 = 0;
601 if ( theShortHitEnergy29 > shortFibre_Cut &&
602 ( theLongHitEnergy29/theShortHitEnergy29 < 2.*longFibre_Fraction ||
607 unsigned theStatusValue = theStatus->
getValue();
609 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
612 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
614 pfrhHFHADCleaned29 = createHcalRecHit( detid,
619 pfrhHFHADCleaned29->
setTime(theShortHit29->time());
633 shortFibre -= theShortHitEnergy29;
634 theShortHitEnergy29 = 0.;
639 if ( theLongHitEnergy29 > longFibre_Cut &&
640 ( theShortHitEnergy29/theLongHitEnergy29 < shortFibre_Fraction ||
644 const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId29);
645 unsigned theStatusValue = theStatus->
getValue();
646 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
649 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
652 pfrhHFEMCleaned29 = createHcalRecHit( detid,
657 pfrhHFEMCleaned29->
setTime(theLongHit29->time());
671 longFibre -= theLongHitEnergy29;
672 theLongHitEnergy29 = 0.;
679 if ( theLongHitEnergy29 >
std::max(theLongHitEnergy,shortFibre_Cut/2) ) {
680 pfrhHFEMCleaned29 = createHcalRecHit( detid,
685 pfrhHFEMCleaned29->
setTime(theLongHit29->time());
693 longFibre -= theLongHitEnergy29;
694 theLongHitEnergy29 = 0.;
697 if ( theShortHitEnergy29 >
std::max(theShortHitEnergy,shortFibre_Cut/2.) ) {
698 pfrhHFHADCleaned29 = createHcalRecHit( detid,
703 pfrhHFHADCleaned29->
setTime(theShortHit29->time());
711 shortFibre -= theShortHitEnergy29;
712 theShortHitEnergy29 = 0.;
718 energyhadHF = 2.*shortFibre;
719 energyemHF = longFibre - shortFibre;
724 if ( energyemHF < thresh_HF_ ) {
725 energyhadHF += energyemHF;
730 if ( HF_Calib_ &&
abs(detid.
ieta()) <= 32 ) {
736 if ( energyemHF > thresh_HF_ || energyhadHF > thresh_HF_ ) {
737 pfrhHFEM = createHcalRecHit( detid,
742 pfrhHFHAD = createHcalRecHit( detid,
754 <<
"CaloTower constituent: unknown layer : " 764 rechitsCleaned->push_back( *pfrhCleaned );
768 HFEMRecHits->push_back( *pfrhHFEM );
772 HFHADRecHits->push_back( *pfrhHFHAD );
775 if(pfrhHFEMCleaned) {
776 rechitsCleaned->push_back( *pfrhHFEMCleaned );
777 delete pfrhHFEMCleaned;
779 if(pfrhHFHADCleaned) {
780 rechitsCleaned->push_back( *pfrhHFHADCleaned );
781 delete pfrhHFHADCleaned;
783 if(pfrhHFEMCleaned29) {
784 rechitsCleaned->push_back( *pfrhHFEMCleaned29 );
785 delete pfrhHFEMCleaned29;
787 if(pfrhHFHADCleaned29) {
788 rechitsCleaned->push_back( *pfrhHFHADCleaned29 );
789 delete pfrhHFHADCleaned29;
801 for(
unsigned int i=0;
i<
rechits->size();++
i) {
805 if (navigation_HF_) {
811 for(
unsigned int i=0;
i<HFEMRecHits->size();++
i) {
812 navigator_->associateNeighbours(HFEMRecHits->at(
i),HFEMRecHits,refProdEM);
819 for(
unsigned int i=0;
i<HFHADRecHits->size();++
i) {
820 navigator_->associateNeighbours(HFHADRecHits->at(
i),HFHADRecHits,refProdHAD);
842 theHcalChStatus = hcalChStatus.
product();
847 theEcalChStatus = ecalChStatus.
product();
851 theTowerConstituentsMap = cttopo.
product();
865 <<
"warning detid "<<detid.
rawId()<<
" not found in layer " 874 auto zp =
dynamic_cast<IdealZPrism const*
>(thisCell);
876 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.
virtual 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 &)
std::vector< CaloTower >::const_iterator const_iterator
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)
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
uint32_t rawId() const
get the raw id
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)
const std::vector< DetId > & constituents() const
void setTime(double time)
RefProd< PROD > getRefBeforePut()
const_iterator end() const
int iphi() const
get the cell iphi
CaloTowerDetId id() const
std::vector< Item >::const_iterator const_iterator
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)