70 shortFibre_Cut = iConfig.
getParameter<
double>(
"ShortFibre_Cut");
71 longFibre_Fraction = iConfig.
getParameter<
double>(
"LongFibre_Fraction");
73 longFibre_Cut = iConfig.
getParameter<
double>(
"LongFibre_Cut");
74 shortFibre_Fraction = iConfig.
getParameter<
double>(
"ShortFibre_Fraction");
76 applyLongShortDPG_ = iConfig.
getParameter<
bool>(
"ApplyLongShortDPG");
78 longShortFibre_Cut = iConfig.
getParameter<
double>(
"LongShortFibre_Cut");
79 minShortTiming_Cut = iConfig.
getParameter<
double>(
"MinShortTiming_Cut");
80 maxShortTiming_Cut = iConfig.
getParameter<
double>(
"MaxShortTiming_Cut");
81 minLongTiming_Cut = iConfig.
getParameter<
double>(
"MinLongTiming_Cut");
82 maxLongTiming_Cut = iConfig.
getParameter<
double>(
"MaxLongTiming_Cut");
84 applyTimeDPG_ = iConfig.
getParameter<
bool>(
"ApplyTimeDPG");
85 applyPulseDPG_ = iConfig.
getParameter<
bool>(
"ApplyPulseDPG");
86 HcalMaxAllowedHFLongShortSev_ = iConfig.
getParameter<
int>(
"HcalMaxAllowedHFLongShortSev");
87 HcalMaxAllowedHFDigiTimeSev_ = iConfig.
getParameter<
int>(
"HcalMaxAllowedHFDigiTimeSev");
88 HcalMaxAllowedHFInTimeWindowSev_ = iConfig.
getParameter<
int>(
"HcalMaxAllowedHFInTimeWindowSev");
89 HcalMaxAllowedChannelStatusSev_ = iConfig.
getParameter<
int>(
"HcalMaxAllowedChannelStatusSev");
91 ECAL_Compensate_ = iConfig.
getParameter<
bool>(
"ECAL_Compensate");
92 ECAL_Threshold_ = iConfig.
getParameter<
double>(
"ECAL_Threshold");
93 ECAL_Compensation_ = iConfig.
getParameter<
double>(
"ECAL_Compensation");
94 ECAL_Dead_Code_ = iConfig.
getParameter<
unsigned int>(
"ECAL_Dead_Code");
113 produces<reco::PFRecHitCollection>();
114 produces<reco::PFRecHitCollection>(
"Cleaned");
115 produces<reco::PFRecHitCollection>(
"HFHAD").setBranchAlias(
"HFHADRecHits");
116 produces<reco::PFRecHitCollection>(
"HFEM").setBranchAlias(
"HFEMRecHits");
125 navigator_->beginEvent(iSetup);
144 auto_ptr< vector<reco::PFRecHit> >
rechits(
new vector<reco::PFRecHit> );
145 auto_ptr< vector<reco::PFRecHit> > rechitsCleaned(
new vector<reco::PFRecHit> );
146 auto_ptr< vector<reco::PFRecHit> > HFHADRecHits(
new vector<reco::PFRecHit> );
147 auto_ptr< vector<reco::PFRecHit> > HFEMRecHits(
new vector<reco::PFRecHit> );
159 for(ICT ict=caloTowers->begin(); ict!=caloTowers->end();ict++) {
173 if( (energy+energyEM) < 1
e-9 )
continue;
177 const std::vector<DetId>& allConstituents = theTowerConstituentsMap->constituentsOf(ct.
id());
179 bool foundHCALConstituent =
false;
183 for(
unsigned int i=0;
i< hits.size();++
i) {
185 foundHCALConstituent =
true;
188 if ( ECAL_Compensate_ && energy > ECAL_Threshold_ ) {
189 for(
unsigned int j=0;
j<allConstituents.size();++
j) {
193 unsigned int dbStatus = chIt != theEcalChStatus->end() ? chIt->getStatusCode() : 0;
194 if ( dbStatus > ECAL_Dead_Code_ ) dead += 1.;
206 double rescaleFactor = alive > 0. ? 1. + ECAL_Compensation_*dead/alive : 1.;
218 if(foundHCALConstituent)
221 switch( detid.
subdet() ) {
224 if(energy < thresh_Barrel_ )
continue;
225 if ( rescaleFactor > 1. ) {
226 pfrhCleaned = createHcalRecHit( detid,
231 pfrhCleaned->
setTime(rescaleFactor);
232 energy *= rescaleFactor;
234 pfrh = createHcalRecHit( detid,
244 if(energy < thresh_Endcap_ )
continue;
247 if ( rescaleFactor > 1. ) {
248 pfrhCleaned = createHcalRecHit( detid,
253 pfrhCleaned->
setTime(rescaleFactor);
254 energy *= rescaleFactor;
256 pfrh = createHcalRecHit( detid,
273 double energyemHF = weight_HFem_ * energyEM;
274 double energyhadHF = weight_HFhad_ *
energy;
276 if((energyemHF+energyhadHF) < thresh_HF_ )
continue;
279 double longFibre = energyemHF + energyhadHF/2.;
280 double shortFibre = energyhadHF/2.;
281 int ieta = detid.
ieta();
282 int iphi = detid.
iphi();
286 iHF theLongHit = hfHandle->find(theLongDetId);
287 iHF theShortHit = hfHandle->find(theShortDetId);
289 double theLongHitEnergy = 0.;
290 double theShortHitEnergy = 0.;
291 bool flagShortDPG =
false;
292 bool flagLongDPG =
false;
293 bool flagShortTimeDPG =
false;
294 bool flagLongTimeDPG =
false;
295 bool flagShortPulseDPG =
false;
296 bool flagLongPulseDPG =
false;
298 if ( theLongHit != hfHandle->end() ) {
299 int theLongFlag = theLongHit->flags();
300 theLongHitEnergy = theLongHit->energy();
301 flagLongDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId, theLongFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
302 flagLongTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId, theLongFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
303 flagLongPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId, theLongFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
310 if ( theShortHit != hfHandle->end() ) {
311 int theShortFlag = theShortHit->flags();
312 theShortHitEnergy = theShortHit->energy();
313 flagShortDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId, theShortFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
314 flagShortTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId, theShortFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
315 flagShortPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId, theShortFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
322 if ( theShortHitEnergy > longShortFibre_Cut &&
323 ( theShortHit->time() < minShortTiming_Cut ||
324 theShortHit->time() > maxShortTiming_Cut ||
325 flagShortTimeDPG || flagShortPulseDPG ) ) {
327 pfrhHFHADCleaned = createHcalRecHit( detid,
332 pfrhHFHADCleaned->
setTime(theShortHit->time());
345 shortFibre -= theShortHitEnergy;
346 theShortHitEnergy = 0.;
350 if ( theLongHitEnergy > longShortFibre_Cut &&
351 ( theLongHit->time() < minLongTiming_Cut ||
352 theLongHit->time() > maxLongTiming_Cut ||
353 flagLongTimeDPG || flagLongPulseDPG ) ) {
355 pfrhHFEMCleaned = createHcalRecHit( detid,
360 pfrhHFEMCleaned->
setTime(theLongHit->time());
373 longFibre -= theLongHitEnergy;
374 theLongHitEnergy = 0.;
378 if ( theShortHitEnergy > shortFibre_Cut &&
379 ( theLongHitEnergy/theShortHitEnergy < longFibre_Fraction ||
384 unsigned theStatusValue = theStatus->
getValue();
385 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
388 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
390 pfrhHFHADCleaned = createHcalRecHit( detid,
395 pfrhHFHADCleaned->
setTime(theShortHit->time());
409 shortFibre -= theShortHitEnergy;
410 theShortHitEnergy = 0.;
414 if ( theLongHitEnergy > longFibre_Cut &&
415 ( theShortHitEnergy/theLongHitEnergy < shortFibre_Fraction ||
420 unsigned theStatusValue = theStatus->
getValue();
422 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
425 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
428 pfrhHFEMCleaned = createHcalRecHit( detid,
433 pfrhHFEMCleaned->
setTime(theLongHit->time());
447 longFibre -= theLongHitEnergy;
448 theLongHitEnergy = 0.;
454 if (
abs(ieta) == 29 ) {
457 if ( theLongHitEnergy > shortFibre_Cut/2. ) {
458 pfrhHFEMCleaned29 = createHcalRecHit( detid,
463 pfrhHFEMCleaned29->
setTime(theLongHit->time());
471 longFibre -= theLongHitEnergy;
472 theLongHitEnergy = 0.;
475 if ( theShortHitEnergy > shortFibre_Cut/2. ) {
476 pfrhHFHADCleaned29 = createHcalRecHit( detid,
481 pfrhHFHADCleaned29->
setTime(theShortHit->time());
489 shortFibre -= theShortHitEnergy;
490 theShortHitEnergy = 0.;
496 else if (
abs(ieta) == 30 ) {
497 int ieta29 = ieta > 0 ? 29 : -29;
500 iHF theLongHit29 = hfHandle->find(theLongDetId29);
501 iHF theShortHit29 = hfHandle->find(theShortDetId29);
503 double theLongHitEnergy29 = 0.;
504 double theShortHitEnergy29 = 0.;
505 bool flagShortDPG29 =
false;
506 bool flagLongDPG29 =
false;
507 bool flagShortTimeDPG29 =
false;
508 bool flagLongTimeDPG29 =
false;
509 bool flagShortPulseDPG29 =
false;
510 bool flagLongPulseDPG29 =
false;
512 if ( theLongHit29 != hfHandle->end() ) {
513 int theLongFlag29 = theLongHit29->flags();
514 theLongHitEnergy29 = theLongHit29->energy() ;
515 flagLongDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
516 flagLongTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
517 flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
524 if ( theShortHit29 != hfHandle->end() ) {
525 int theShortFlag29 = theShortHit29->flags();
526 theShortHitEnergy29 = theShortHit29->energy();
527 flagShortDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
528 flagShortTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
529 flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
536 if ( theLongHitEnergy29 > longShortFibre_Cut &&
537 ( theLongHit29->time() < minLongTiming_Cut ||
538 theLongHit29->time() > maxLongTiming_Cut ||
539 flagLongTimeDPG29 || flagLongPulseDPG29 ) ) {
541 pfrhHFEMCleaned29 = createHcalRecHit( detid,
546 pfrhHFEMCleaned29->
setTime(theLongHit29->time());
559 longFibre -= theLongHitEnergy29;
560 theLongHitEnergy29 = 0;
563 if ( theShortHitEnergy29 > longShortFibre_Cut &&
564 ( theShortHit29->time() < minShortTiming_Cut ||
565 theShortHit29->time() > maxShortTiming_Cut ||
566 flagShortTimeDPG29 || flagShortPulseDPG29 ) ) {
568 pfrhHFHADCleaned29 = createHcalRecHit( detid,
573 pfrhHFHADCleaned29->
setTime(theShortHit29->time());
586 shortFibre -= theShortHitEnergy29;
587 theShortHitEnergy29 = 0.;
591 if ( theShortHitEnergy29 > shortFibre_Cut &&
592 ( theLongHitEnergy29/theShortHitEnergy29 < 2.*longFibre_Fraction ||
597 unsigned theStatusValue = theStatus->
getValue();
599 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
602 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
604 pfrhHFHADCleaned29 = createHcalRecHit( detid,
609 pfrhHFHADCleaned29->
setTime(theShortHit29->time());
623 shortFibre -= theShortHitEnergy29;
624 theShortHitEnergy29 = 0.;
629 if ( theLongHitEnergy29 > longFibre_Cut &&
630 ( theShortHitEnergy29/theLongHitEnergy29 < shortFibre_Fraction ||
634 const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId29);
635 unsigned theStatusValue = theStatus->
getValue();
636 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
639 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
642 pfrhHFEMCleaned29 = createHcalRecHit( detid,
647 pfrhHFEMCleaned29->
setTime(theLongHit29->time());
661 longFibre -= theLongHitEnergy29;
662 theLongHitEnergy29 = 0.;
669 if ( theLongHitEnergy29 >
std::max(theLongHitEnergy,shortFibre_Cut/2) ) {
670 pfrhHFEMCleaned29 = createHcalRecHit( detid,
675 pfrhHFEMCleaned29->
setTime(theLongHit29->time());
683 longFibre -= theLongHitEnergy29;
684 theLongHitEnergy29 = 0.;
687 if ( theShortHitEnergy29 >
std::max(theShortHitEnergy,shortFibre_Cut/2.) ) {
688 pfrhHFHADCleaned29 = createHcalRecHit( detid,
693 pfrhHFHADCleaned29->
setTime(theShortHit29->time());
701 shortFibre -= theShortHitEnergy29;
702 theShortHitEnergy29 = 0.;
708 energyhadHF = 2.*shortFibre;
709 energyemHF = longFibre - shortFibre;
714 if ( energyemHF < thresh_HF_ ) {
715 energyhadHF += energyemHF;
720 if ( HF_Calib_ &&
abs(detid.
ieta()) <= 32 ) {
726 if ( energyemHF > thresh_HF_ || energyhadHF > thresh_HF_ ) {
727 pfrhHFEM = createHcalRecHit( detid,
732 pfrhHFHAD = createHcalRecHit( detid,
744 <<
"CaloTower constituent: unknown layer : "
750 rechits->push_back( *pfrh );
754 rechitsCleaned->push_back( *pfrhCleaned );
758 HFEMRecHits->push_back( *pfrhHFEM );
762 HFHADRecHits->push_back( *pfrhHFHAD );
765 if(pfrhHFEMCleaned) {
766 rechitsCleaned->push_back( *pfrhHFEMCleaned );
767 delete pfrhHFEMCleaned;
769 if(pfrhHFHADCleaned) {
770 rechitsCleaned->push_back( *pfrhHFHADCleaned );
771 delete pfrhHFHADCleaned;
773 if(pfrhHFEMCleaned29) {
774 rechitsCleaned->push_back( *pfrhHFEMCleaned29 );
775 delete pfrhHFEMCleaned29;
777 if(pfrhHFHADCleaned29) {
778 rechitsCleaned->push_back( *pfrhHFHADCleaned29 );
779 delete pfrhHFHADCleaned29;
791 for(
unsigned int i=0;
i<rechits->size();++
i) {
792 navigator_->associateNeighbours(rechits->at(
i),
rechits,refProd);
795 if (navigation_HF_) {
801 for(
unsigned int i=0;
i<HFEMRecHits->size();++
i) {
802 navigator_->associateNeighbours(HFEMRecHits->at(
i),HFEMRecHits,refProdEM);
809 for(
unsigned int i=0;
i<HFHADRecHits->size();++
i) {
810 navigator_->associateNeighbours(HFHADRecHits->at(
i),HFHADRecHits,refProdHAD);
814 iEvent.
put( rechits,
"" );
815 iEvent.
put( rechitsCleaned,
"Cleaned" );
816 iEvent.
put( HFEMRecHits,
"HFEM" );
817 iEvent.
put( HFHADRecHits,
"HFHAD" );
832 theHcalChStatus = hcalChStatus.
product();
837 theEcalChStatus = ecalChStatus.
product();
841 theTowerConstituentsMap = cttopo.
product();
855 <<
"warning detid "<<detid.
rawId()<<
" not found in layer "
862 double depth_correction = 0.;
865 depth_correction = position.
z() > 0. ? EM_Depth_ : -EM_Depth_;
868 depth_correction = position.
z() > 0. ? HAD_Depth_ : -HAD_Depth_;
876 position.
x(), position.
y(), position.
z()+depth_correction,
885 rh->
setNECorner( corners[0].
x(), corners[0].y(), corners[0].z()+depth_correction );
886 rh->
setSECorner( corners[1].
x(), corners[1].y(), corners[1].z()+depth_correction );
887 rh->
setSWCorner( corners[2].
x(), corners[2].y(), corners[2].z()+depth_correction );
888 rh->
setNWCorner( corners[3].
x(), corners[3].y(), corners[3].z()+depth_correction );
void setSECorner(double posx, double posy, double posz)
T getParameter(std::string const &) const
void setNECorner(double posx, double posy, double posz)
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
void setSWCorner(double posx, double posy, double posz)
T x() const
Cartesian x coordinate.
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
void produce(edm::Event &iEvent, const edm::EventSetup &iSetup)
int ieta() const
get the cell ieta
Abs< T >::type abs(const T &t)
const std::vector< DetId > & constituents() const
void setTime(double time)
void setNWCorner(double posx, double posy, double posz)
RefProd< PROD > getRefBeforePut()
int iphi() const
get the cell iphi
CaloTowerDetId id() const
std::vector< Item >::const_iterator const_iterator
T const * product() const
int getSeverityLevel(const DetId &myid, const uint32_t &myflag, const uint32_t &mystatus) const
static int position[264][3]
uint32_t getValue() const
const CornersVec & getCorners() const
Returns the corner points of this cell's volume.
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
T get(const Candidate &c)