71 shortFibre_Cut = iConfig.
getParameter<
double>(
"ShortFibre_Cut");
72 longFibre_Fraction = iConfig.
getParameter<
double>(
"LongFibre_Fraction");
74 longFibre_Cut = iConfig.
getParameter<
double>(
"LongFibre_Cut");
75 shortFibre_Fraction = iConfig.
getParameter<
double>(
"ShortFibre_Fraction");
77 applyLongShortDPG_ = iConfig.
getParameter<
bool>(
"ApplyLongShortDPG");
79 longShortFibre_Cut = iConfig.
getParameter<
double>(
"LongShortFibre_Cut");
80 minShortTiming_Cut = iConfig.
getParameter<
double>(
"MinShortTiming_Cut");
81 maxShortTiming_Cut = iConfig.
getParameter<
double>(
"MaxShortTiming_Cut");
82 minLongTiming_Cut = iConfig.
getParameter<
double>(
"MinLongTiming_Cut");
83 maxLongTiming_Cut = iConfig.
getParameter<
double>(
"MaxLongTiming_Cut");
85 applyTimeDPG_ = iConfig.
getParameter<
bool>(
"ApplyTimeDPG");
86 applyPulseDPG_ = iConfig.
getParameter<
bool>(
"ApplyPulseDPG");
87 HcalMaxAllowedHFLongShortSev_ = iConfig.
getParameter<
int>(
"HcalMaxAllowedHFLongShortSev");
88 HcalMaxAllowedHFDigiTimeSev_ = iConfig.
getParameter<
int>(
"HcalMaxAllowedHFDigiTimeSev");
89 HcalMaxAllowedHFInTimeWindowSev_ = iConfig.
getParameter<
int>(
"HcalMaxAllowedHFInTimeWindowSev");
90 HcalMaxAllowedChannelStatusSev_ = iConfig.
getParameter<
int>(
"HcalMaxAllowedChannelStatusSev");
92 ECAL_Compensate_ = iConfig.
getParameter<
bool>(
"ECAL_Compensate");
93 ECAL_Threshold_ = iConfig.
getParameter<
double>(
"ECAL_Threshold");
94 ECAL_Compensation_ = iConfig.
getParameter<
double>(
"ECAL_Compensation");
95 ECAL_Dead_Code_ = iConfig.
getParameter<
unsigned int>(
"ECAL_Dead_Code");
114 produces<reco::PFRecHitCollection>();
115 produces<reco::PFRecHitCollection>(
"Cleaned");
116 produces<reco::PFRecHitCollection>(
"HFHAD").setBranchAlias(
"HFHADRecHits");
117 produces<reco::PFRecHitCollection>(
"HFEM").setBranchAlias(
"HFEMRecHits");
126 navigator_->beginEvent(iSetup);
145 auto_ptr< vector<reco::PFRecHit> >
rechits(
new vector<reco::PFRecHit> );
146 auto_ptr< vector<reco::PFRecHit> > rechitsCleaned(
new vector<reco::PFRecHit> );
147 auto_ptr< vector<reco::PFRecHit> > HFHADRecHits(
new vector<reco::PFRecHit> );
148 auto_ptr< vector<reco::PFRecHit> > HFEMRecHits(
new vector<reco::PFRecHit> );
160 for(ICT ict=caloTowers->begin(); ict!=caloTowers->end();ict++) {
174 if( (energy+energyEM) < 1
e-9 )
continue;
178 const std::vector<DetId>& allConstituents = theTowerConstituentsMap->constituentsOf(ct.
id());
180 bool foundHCALConstituent =
false;
184 for(
unsigned int i=0;
i< hits.size();++
i) {
186 foundHCALConstituent =
true;
189 if ( ECAL_Compensate_ && energy > ECAL_Threshold_ ) {
190 for(
unsigned int j=0;
j<allConstituents.size();++
j) {
194 unsigned int dbStatus = chIt != theEcalChStatus->end() ? chIt->getStatusCode() : 0;
195 if ( dbStatus > ECAL_Dead_Code_ ) dead += 1.;
207 double rescaleFactor = alive > 0. ? 1. + ECAL_Compensation_*dead/alive : 1.;
219 if(foundHCALConstituent)
222 switch( detid.
subdet() ) {
225 if(energy < thresh_Barrel_ )
continue;
226 if ( rescaleFactor > 1. ) {
227 pfrhCleaned = createHcalRecHit( detid,
232 pfrhCleaned->
setTime(rescaleFactor);
233 energy *= rescaleFactor;
235 pfrh = createHcalRecHit( detid,
245 if(energy < thresh_Endcap_ )
continue;
247 if ( HCAL_Calib_ &&
abs(detid.
ieta()) == 29 ) energy *= HCAL_Calib_29;
248 if ( rescaleFactor > 1. ) {
249 pfrhCleaned = createHcalRecHit( detid,
254 pfrhCleaned->
setTime(rescaleFactor);
255 energy *= rescaleFactor;
257 pfrh = createHcalRecHit( detid,
274 double energyemHF = weight_HFem_ * energyEM;
275 double energyhadHF = weight_HFhad_ *
energy;
277 if((energyemHF+energyhadHF) < thresh_HF_ )
continue;
280 double longFibre = energyemHF + energyhadHF/2.;
281 double shortFibre = energyhadHF/2.;
282 int ieta = detid.
ieta();
283 int iphi = detid.
iphi();
287 iHF theLongHit = hfHandle->find(theLongDetId);
288 iHF theShortHit = hfHandle->find(theShortDetId);
290 double theLongHitEnergy = 0.;
291 double theShortHitEnergy = 0.;
292 bool flagShortDPG =
false;
293 bool flagLongDPG =
false;
294 bool flagShortTimeDPG =
false;
295 bool flagLongTimeDPG =
false;
296 bool flagShortPulseDPG =
false;
297 bool flagLongPulseDPG =
false;
299 if ( theLongHit != hfHandle->end() ) {
300 int theLongFlag = theLongHit->flags();
301 theLongHitEnergy = theLongHit->energy();
302 flagLongDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId, theLongFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
303 flagLongTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId, theLongFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
304 flagLongPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId, theLongFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
311 if ( theShortHit != hfHandle->end() ) {
312 int theShortFlag = theShortHit->flags();
313 theShortHitEnergy = theShortHit->energy();
314 flagShortDPG = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId, theShortFlag & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
315 flagShortTimeDPG = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId, theShortFlag & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
316 flagShortPulseDPG = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId, theShortFlag & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
323 if ( theShortHitEnergy > longShortFibre_Cut &&
324 ( theShortHit->time() < minShortTiming_Cut ||
325 theShortHit->time() > maxShortTiming_Cut ||
326 flagShortTimeDPG || flagShortPulseDPG ) ) {
328 pfrhHFHADCleaned = createHcalRecHit( detid,
333 pfrhHFHADCleaned->
setTime(theShortHit->time());
346 shortFibre -= theShortHitEnergy;
347 theShortHitEnergy = 0.;
351 if ( theLongHitEnergy > longShortFibre_Cut &&
352 ( theLongHit->time() < minLongTiming_Cut ||
353 theLongHit->time() > maxLongTiming_Cut ||
354 flagLongTimeDPG || flagLongPulseDPG ) ) {
356 pfrhHFEMCleaned = createHcalRecHit( detid,
361 pfrhHFEMCleaned->
setTime(theLongHit->time());
374 longFibre -= theLongHitEnergy;
375 theLongHitEnergy = 0.;
379 if ( theShortHitEnergy > shortFibre_Cut &&
380 ( theLongHitEnergy/theShortHitEnergy < longFibre_Fraction ||
385 unsigned theStatusValue = theStatus->
getValue();
386 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
389 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
391 pfrhHFHADCleaned = createHcalRecHit( detid,
396 pfrhHFHADCleaned->
setTime(theShortHit->time());
410 shortFibre -= theShortHitEnergy;
411 theShortHitEnergy = 0.;
415 if ( theLongHitEnergy > longFibre_Cut &&
416 ( theShortHitEnergy/theLongHitEnergy < shortFibre_Fraction ||
421 unsigned theStatusValue = theStatus->
getValue();
423 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
426 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
429 pfrhHFEMCleaned = createHcalRecHit( detid,
434 pfrhHFEMCleaned->
setTime(theLongHit->time());
448 longFibre -= theLongHitEnergy;
449 theLongHitEnergy = 0.;
455 if (
abs(ieta) == 29 ) {
458 if ( theLongHitEnergy > shortFibre_Cut/2. ) {
459 pfrhHFEMCleaned29 = createHcalRecHit( detid,
464 pfrhHFEMCleaned29->
setTime(theLongHit->time());
472 longFibre -= theLongHitEnergy;
473 theLongHitEnergy = 0.;
476 if ( theShortHitEnergy > shortFibre_Cut/2. ) {
477 pfrhHFHADCleaned29 = createHcalRecHit( detid,
482 pfrhHFHADCleaned29->
setTime(theShortHit->time());
490 shortFibre -= theShortHitEnergy;
491 theShortHitEnergy = 0.;
497 else if (
abs(ieta) == 30 ) {
498 int ieta29 = ieta > 0 ? 29 : -29;
501 iHF theLongHit29 = hfHandle->find(theLongDetId29);
502 iHF theShortHit29 = hfHandle->find(theShortDetId29);
504 double theLongHitEnergy29 = 0.;
505 double theShortHitEnergy29 = 0.;
506 bool flagShortDPG29 =
false;
507 bool flagLongDPG29 =
false;
508 bool flagShortTimeDPG29 =
false;
509 bool flagLongTimeDPG29 =
false;
510 bool flagShortPulseDPG29 =
false;
511 bool flagLongPulseDPG29 =
false;
513 if ( theLongHit29 != hfHandle->end() ) {
514 int theLongFlag29 = theLongHit29->flags();
515 theLongHitEnergy29 = theLongHit29->energy() ;
516 flagLongDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
517 flagLongTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
518 flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theLongDetId29, theLongFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
525 if ( theShortHit29 != hfHandle->end() ) {
526 int theShortFlag29 = theShortHit29->flags();
527 theShortHitEnergy29 = theShortHit29->energy();
528 flagShortDPG29 = applyLongShortDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFLongShortFlagValue_, 0)> HcalMaxAllowedHFLongShortSev_);
529 flagShortTimeDPG29 = applyTimeDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFInTimeWindowFlagValue_, 0)> HcalMaxAllowedHFInTimeWindowSev_);
530 flagLongPulseDPG29 = applyPulseDPG_ && ( hcalSevLvlComputer->
getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFDigiTimeFlagValue_, 0)> HcalMaxAllowedHFDigiTimeSev_);
537 if ( theLongHitEnergy29 > longShortFibre_Cut &&
538 ( theLongHit29->time() < minLongTiming_Cut ||
539 theLongHit29->time() > maxLongTiming_Cut ||
540 flagLongTimeDPG29 || flagLongPulseDPG29 ) ) {
542 pfrhHFEMCleaned29 = createHcalRecHit( detid,
547 pfrhHFEMCleaned29->
setTime(theLongHit29->time());
560 longFibre -= theLongHitEnergy29;
561 theLongHitEnergy29 = 0;
564 if ( theShortHitEnergy29 > longShortFibre_Cut &&
565 ( theShortHit29->time() < minShortTiming_Cut ||
566 theShortHit29->time() > maxShortTiming_Cut ||
567 flagShortTimeDPG29 || flagShortPulseDPG29 ) ) {
569 pfrhHFHADCleaned29 = createHcalRecHit( detid,
574 pfrhHFHADCleaned29->
setTime(theShortHit29->time());
587 shortFibre -= theShortHitEnergy29;
588 theShortHitEnergy29 = 0.;
592 if ( theShortHitEnergy29 > shortFibre_Cut &&
593 ( theLongHitEnergy29/theShortHitEnergy29 < 2.*longFibre_Fraction ||
598 unsigned theStatusValue = theStatus->
getValue();
600 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
603 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
605 pfrhHFHADCleaned29 = createHcalRecHit( detid,
610 pfrhHFHADCleaned29->
setTime(theShortHit29->time());
624 shortFibre -= theShortHitEnergy29;
625 theShortHitEnergy29 = 0.;
630 if ( theLongHitEnergy29 > longFibre_Cut &&
631 ( theShortHitEnergy29/theLongHitEnergy29 < shortFibre_Fraction ||
635 const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId29);
636 unsigned theStatusValue = theStatus->
getValue();
637 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
640 if (theSeverityLevel<=HcalMaxAllowedChannelStatusSev_) {
643 pfrhHFEMCleaned29 = createHcalRecHit( detid,
648 pfrhHFEMCleaned29->
setTime(theLongHit29->time());
662 longFibre -= theLongHitEnergy29;
663 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 ) {
721 energyhadHF *= HF_Calib_29;
722 energyemHF *= HF_Calib_29;
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)
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)