115 produces<reco::PFRecHitCollection>(
"HFHAD").setBranchAlias(
"HFHADRecHits");
116 produces<reco::PFRecHitCollection>(
"HFEM").setBranchAlias(
"HFEMRecHits");
127 vector<reco::PFRecHit>& rechitsCleaned,
139 map<unsigned, unsigned > idSortedRecHits;
140 map<unsigned, unsigned > idSortedRecHitsHFEM;
141 map<unsigned, unsigned > idSortedRecHitsHFHAD;
142 typedef map<unsigned, unsigned >::iterator
IDH;
162 auto_ptr< vector<reco::PFRecHit> > HFHADRecHits(
new vector<reco::PFRecHit> );
163 auto_ptr< vector<reco::PFRecHit> > HFEMRecHits(
new vector<reco::PFRecHit> );
192 LogError(
"PFRecHitProducerHCAL")<<err.str()<<endl;
197 assert( caloTowers.
isValid() );
207 LogError(
"PFRecHitProducerHCAL")<<err.str()<<endl;
212 assert( hfHandle.isValid() );
223 LogError(
"PFRecHitProducerHCAL")<<err.str()<<endl;
228 assert( hbheHandle.
isValid() );
234 for(ICT ict=caloTowers->begin(); ict!=caloTowers->end();ict++) {
239 if(!caloTowerGeometry)
240 caloTowerGeometry = geoHandle->getSubdetectorGeometry(ct.
id());
253 if( (energy+energyEM) < 1
e-9 )
continue;
291 bool foundHCALConstituent =
false;
297 for(
unsigned int i=0;
i< hits.size();++
i) {
299 foundHCALConstituent =
true;
303 for(
unsigned int j=0;
j<allConstituents.size();++
j) {
332 if(foundHCALConstituent)
335 switch( detid.
subdet() ) {
360 if ( rescaleFactor > 1. ) {
367 energy *= rescaleFactor;
401 if ( rescaleFactor > 1. ) {
408 energy *= rescaleFactor;
430 if((energyemHF+energyhadHF) <
thresh_HF_ )
continue;
433 double longFibre = energyemHF + energyhadHF/2.;
434 double shortFibre = energyhadHF/2.;
435 int ieta = detid.
ieta();
436 int iphi = detid.
iphi();
440 iHF theLongHit = hfHandle->find(theLongDetId);
441 iHF theShortHit = hfHandle->find(theShortDetId);
443 double theLongHitEnergy = 0.;
444 double theShortHitEnergy = 0.;
445 bool flagShortDPG =
false;
446 bool flagLongDPG =
false;
447 bool flagShortTimeDPG =
false;
448 bool flagLongTimeDPG =
false;
449 bool flagShortPulseDPG =
false;
450 bool flagLongPulseDPG =
false;
452 if ( theLongHit != hfHandle->end() ) {
453 int theLongFlag = theLongHit->flags();
454 theLongHitEnergy = theLongHit->energy();
464 if ( theShortHit != hfHandle->end() ) {
465 int theShortFlag = theShortHit->flags();
466 theShortHitEnergy = theShortHit->energy();
479 flagShortTimeDPG || flagShortPulseDPG ) ) {
486 pfrhHFHADCleaned->
setRescale(theShortHit->time());
499 shortFibre -= theShortHitEnergy;
500 theShortHitEnergy = 0.;
507 flagLongTimeDPG || flagLongPulseDPG ) ) {
514 pfrhHFEMCleaned->
setRescale(theLongHit->time());
527 longFibre -= theLongHitEnergy;
528 theLongHitEnergy = 0.;
538 unsigned theStatusValue = theStatus->
getValue();
539 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
549 pfrhHFHADCleaned->
setRescale(theShortHit->time());
563 shortFibre -= theShortHitEnergy;
564 theShortHitEnergy = 0.;
574 unsigned theStatusValue = theStatus->
getValue();
576 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
587 pfrhHFEMCleaned->
setRescale(theLongHit->time());
601 longFibre -= theLongHitEnergy;
602 theLongHitEnergy = 0.;
608 if (
abs(ieta) == 29 ) {
617 pfrhHFEMCleaned29->
setRescale(theLongHit->time());
625 longFibre -= theLongHitEnergy;
626 theLongHitEnergy = 0.;
635 pfrhHFHADCleaned29->
setRescale(theShortHit->time());
643 shortFibre -= theShortHitEnergy;
644 theShortHitEnergy = 0.;
650 else if (
abs(ieta) == 30 ) {
651 int ieta29 = ieta > 0 ? 29 : -29;
654 iHF theLongHit29 = hfHandle->find(theLongDetId29);
655 iHF theShortHit29 = hfHandle->find(theShortDetId29);
657 double theLongHitEnergy29 = 0.;
658 double theShortHitEnergy29 = 0.;
659 bool flagShortDPG29 =
false;
660 bool flagLongDPG29 =
false;
661 bool flagShortTimeDPG29 =
false;
662 bool flagLongTimeDPG29 =
false;
663 bool flagShortPulseDPG29 =
false;
664 bool flagLongPulseDPG29 =
false;
666 if ( theLongHit29 != hfHandle->end() ) {
667 int theLongFlag29 = theLongHit29->flags();
668 theLongHitEnergy29 = theLongHit29->energy() ;
678 if ( theShortHit29 != hfHandle->end() ) {
679 int theShortFlag29 = theShortHit29->flags();
680 theShortHitEnergy29 = theShortHit29->energy();
693 flagLongTimeDPG29 || flagLongPulseDPG29 ) ) {
700 pfrhHFEMCleaned29->
setRescale(theLongHit29->time());
713 longFibre -= theLongHitEnergy29;
714 theLongHitEnergy29 = 0;
720 flagShortTimeDPG29 || flagShortPulseDPG29 ) ) {
727 pfrhHFHADCleaned29->
setRescale(theShortHit29->time());
740 shortFibre -= theShortHitEnergy29;
741 theShortHitEnergy29 = 0.;
751 unsigned theStatusValue = theStatus->
getValue();
753 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
763 pfrhHFHADCleaned29->
setRescale(theShortHit29->time());
777 shortFibre -= theShortHitEnergy29;
778 theShortHitEnergy29 = 0.;
789 unsigned theStatusValue = theStatus->
getValue();
790 int theSeverityLevel = hcalSevLvlComputer->
getSeverityLevel(detid, 0, theStatusValue);
801 pfrhHFEMCleaned29->
setRescale(theLongHit29->time());
815 longFibre -= theLongHitEnergy29;
816 theLongHitEnergy29 = 0.;
828 pfrhHFEMCleaned29->
setRescale(theLongHit29->time());
836 longFibre -= theLongHitEnergy29;
837 theLongHitEnergy29 = 0.;
846 pfrhHFHADCleaned29->
setRescale(theShortHit29->time());
854 shortFibre -= theShortHitEnergy29;
855 theShortHitEnergy29 = 0.;
861 energyhadHF = 2.*shortFibre;
862 energyemHF = longFibre - shortFibre;
868 energyhadHF += energyemHF;
891 pfrhHFHAD->setEnergyUp(energyemHF);
898 <<
"CaloTower constituent: unknown layer : "
903 rechits.push_back( *pfrh );
905 idSortedRecHits.insert( make_pair(ct.
id().
rawId(),
906 rechits.size()-1 ) );
909 rechitsCleaned.push_back( *pfrhCleaned );
914 HFEMRecHits->push_back( *pfrhHFEM );
916 idSortedRecHitsHFEM.insert( make_pair(ct.
id().
rawId(),
917 HFEMRecHits->size()-1 ) );
920 HFHADRecHits->push_back( *pfrhHFHAD );
922 idSortedRecHitsHFHAD.insert( make_pair(ct.
id().
rawId(),
923 HFHADRecHits->size()-1 ) );
926 if(pfrhHFEMCleaned) {
927 rechitsCleaned.push_back( *pfrhHFEMCleaned );
928 delete pfrhHFEMCleaned;
930 if(pfrhHFHADCleaned) {
931 rechitsCleaned.push_back( *pfrhHFHADCleaned );
932 delete pfrhHFHADCleaned;
934 if(pfrhHFEMCleaned29) {
935 rechitsCleaned.push_back( *pfrhHFEMCleaned29 );
936 delete pfrhHFEMCleaned29;
938 if(pfrhHFHADCleaned29) {
939 rechitsCleaned.push_back( *pfrhHFHADCleaned29 );
940 delete pfrhHFHADCleaned29;
945 for(
unsigned i=0;
i<rechits.size();
i++ ) {
950 for(
unsigned i=0;
i<HFEMRecHits->size();
i++ ) {
955 for(
unsigned i=0;
i<HFHADRecHits->size();
i++ ) {
957 idSortedRecHitsHFHAD,
960 iEvent.
put( HFHADRecHits,
"HFHAD" );
961 iEvent.
put( HFEMRecHits,
"HFEM" );
982 LogError(
"PFRecHitProducerHCAL")<<err.str()<<endl;
987 assert( hcalHandle.
isValid() );
990 for(
unsigned irechit=0; irechit<handle->size(); irechit++) {
999 switch( detid.
subdet() ) {
1006 hcalBarrelGeometry );
1015 hcalEndcapGeometry );
1024 hcalEndcapGeometry );
1029 <<
"HCAL rechit: unknown layer : "<<detid.
subdet()<<endl;
1034 rechits.push_back( *pfrh );
1036 idSortedRecHits.insert( make_pair(detid.
rawId(),
1037 rechits.size()-1 ) );
1043 for(
unsigned i=0;
i<rechits.size();
i++ ) {
1047 *hcalBarrelGeometry,
1049 *hcalEndcapGeometry);
1065 unsigned newDetId ) {
1070 <<
"warning detid "<<detid.
rawId()<<
" not found in layer "
1077 double depth_correction = 0.;
1089 unsigned id =
detid;
1090 if(newDetId)
id = newDetId;
1093 position.
x(), position.
y(), position.
z()+depth_correction,
1102 assert( corners.size() == 8 );
1104 rh->
setNECorner( corners[0].
x(), corners[0].
y(), corners[0].
z()+depth_correction );
1105 rh->
setSECorner( corners[1].
x(), corners[1].
y(), corners[1].
z()+depth_correction );
1106 rh->
setSWCorner( corners[2].
x(), corners[2].
y(), corners[2].
z()+depth_correction );
1107 rh->
setNWCorner( corners[3].
x(), corners[3].
y(), corners[3].
z()+depth_correction );
1118 const map<unsigned,unsigned >& sortedHits,
1125 if(navigation_HF_ ==
false){
1136 switch( rh.
layer() ) {
1138 topology = &endcapTopology;
1139 geometry = &endcapGeometry;
1142 topology = &barrelTopology;
1143 geometry = &barrelGeometry;
1146 topology = &endcapTopology;
1147 geometry = &endcapGeometry;
1150 topology = &barrelTopology;
1151 geometry = &barrelGeometry;
1155 topology = &barrelTopology;
1156 geometry = &barrelGeometry;
1162 assert( topology && geometry );
1169 if( north !=
DetId(0) ) {
1170 northeast = navigator.
east();
1180 if( south !=
DetId(0) ) {
1181 southwest = navigator.
west();
1188 if( east !=
DetId(0) ) {
1189 southeast = navigator.
south();
1194 if( west !=
DetId(0) ) {
1195 northwest = navigator.
north();
1199 IDH i = sortedHits.find( north.
rawId() );
1200 if(i != sortedHits.end() )
1203 i = sortedHits.find( northeast.
rawId() );
1204 if(i != sortedHits.end() )
1207 i = sortedHits.find( south.
rawId() );
1208 if(i != sortedHits.end() )
1211 i = sortedHits.find( southwest.
rawId() );
1212 if(i != sortedHits.end() )
1215 i = sortedHits.find( east.
rawId() );
1216 if(i != sortedHits.end() )
1219 i = sortedHits.find( southeast.
rawId() );
1220 if(i != sortedHits.end() )
1223 i = sortedHits.find( west.
rawId() );
1224 if(i != sortedHits.end() )
1227 i = sortedHits.find( northwest.
rawId() );
1228 if(i != sortedHits.end() )
1238 const map<unsigned, unsigned >& sortedHits,
1246 if(navigation_HF_ ==
false){
1255 vector<DetId> northids = topology.
north(ctDetId);
1256 vector<DetId> westids = topology.
west(ctDetId);
1257 vector<DetId> southids = topology.
south(ctDetId);
1258 vector<DetId> eastids = topology.
east(ctDetId);
1281 switch( northids.size() ) {
1285 north = northids[0];
1288 stringstream err(
"PFRecHitProducerHCAL::findRecHitNeighboursCT : incorrect number of neighbours north: ");
1289 err<<northids.size();
1293 switch( southids.size() ) {
1297 south = southids[0];
1300 stringstream err(
"PFRecHitProducerHCAL::findRecHitNeighboursCT : incorrect number of neighbours south: ");
1301 err<<southids.size();
1308 switch( eastids.size() ) {
1313 northeast = getNorth(east, topology);
1314 southeast = getSouth(east, topology);
1320 northeast = getNorth(east, topology );
1321 southeast = getSouth(east2, topology);
1322 northeast2 = getNorth(northeast, topology );
1323 southeast2 = getSouth(southeast, topology);
1326 stringstream err(
"PFRecHitProducerHCAL::findRecHitNeighboursCT : incorrect number of neighbours eastids: ");
1327 err<<eastids.size();
1332 switch( westids.size() ) {
1337 northwest = getNorth(west, topology);
1338 southwest = getSouth(west, topology);
1344 northwest = getNorth(west, topology );
1345 southwest = getSouth(west2, topology );
1346 northwest2 = getNorth(northwest, topology );
1347 southwest2 = getSouth(southwest, topology );
1350 stringstream err(
"PFRecHitProducerHCAL::findRecHitNeighboursCT : incorrect number of neighbours westids: ");
1351 err<< westids.size();
1360 IDH i = sortedHits.find( north.
rawId() );
1361 if(i != sortedHits.end() )
1364 i = sortedHits.find( northeast.
rawId() );
1365 if(i != sortedHits.end() )
1368 i = sortedHits.find( northeast2.
rawId() );
1369 if(i != sortedHits.end() )
1372 i = sortedHits.find( south.
rawId() );
1373 if(i != sortedHits.end() )
1376 i = sortedHits.find( southwest.
rawId() );
1377 if(i != sortedHits.end() )
1380 i = sortedHits.find( southwest2.
rawId() );
1381 if(i != sortedHits.end() )
1384 i = sortedHits.find( east.
rawId() );
1385 if(i != sortedHits.end() )
1388 i = sortedHits.find( east2.
rawId() );
1389 if(i != sortedHits.end() )
1392 i = sortedHits.find( southeast.
rawId() );
1393 if(i != sortedHits.end() )
1396 i = sortedHits.find( southeast2.
rawId() );
1397 if(i != sortedHits.end() )
1400 i = sortedHits.find( west.
rawId() );
1401 if(i != sortedHits.end() )
1404 i = sortedHits.find( west2.
rawId() );
1405 if(i != sortedHits.end() )
1408 i = sortedHits.find( northwest.
rawId() );
1409 if(i != sortedHits.end() )
1412 i = sortedHits.find( northwest2.
rawId() );
1413 if(i != sortedHits.end() )
1430 vector<DetId> sids = topology.
south(
id);
1431 if(sids.size() == 1)
1444 vector<DetId> nids = topology.
north(
id);
1445 if(nids.size() == 1)
void setSECorner(double posx, double posy, double posz)
T getParameter(std::string const &) const
void add4Neighbour(unsigned index)
int HcalMaxAllowedHFDigiTimeSev_
void add8Neighbour(unsigned index)
size_t constituentsSize() const
void setNECorner(double posx, double posy, double posz)
HcalSubdetector subdet() const
get the subdetector
void home() const
move the navigator back to the starting point
const DetId & detid() const
unsigned detId() const
rechit detId
int HcalMaxAllowedChannelStatusSev_
unsigned int ECAL_Dead_Code_
virtual T west() const
move the navigator west
std::vector< CaloTower >::const_iterator const_iterator
virtual std::vector< DetId > west(const DetId &id) const =0
int hcalHFLongShortFlagValue_
reco::PFRecHit * createHcalRecHit(const DetId &detid, double energy, PFLayer::Layer layer, const CaloSubdetectorGeometry *geom, unsigned newDetId=0)
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
const EcalChannelStatus * theEcalChStatus
int HcalMaxAllowedHFLongShortSev_
virtual T north() const
move the navigator north
static int position[TOTALCHAMBERS][3]
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
std::vector< DetId > constituentsOf(const CaloTowerDetId &id) const
Get the constituent detids for this tower id ( not yet implemented )
double minShortTiming_Cut
double ECAL_Compensation_
uint32_t rawId() const
get the raw id
PFLayer::Layer layer() const
rechit layer
void setSWCorner(double posx, double posy, double posz)
double longShortFibre_Cut
edm::InputTag inputTagCaloTowers_
virtual T east() const
move the navigator east
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Base producer for particle flow rechits (PFRecHit)
const T & max(const T &a, const T &b)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
double shortFibre_Fraction
void setEnergyUp(double eUp)
edm::InputTag inputTagHcalRecHitsHF_
int ieta() const
get the cell ieta
const std::vector< DetId > & constituents() const
DetId getSouth(const DetId &id, const CaloSubdetectorTopology &topology)
void setRescale(double factor)
double longFibre_Fraction
void setNWCorner(double posx, double posy, double posz)
search for pointers to neighbours, using neighbours' DetId.
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
edm::InputTag inputTagHcalRecHitsHBHE_
void findRecHitNeighbours(reco::PFRecHit &rh, const std::map< unsigned, unsigned > &sortedHits, const CaloSubdetectorTopology &barrelTopo, const CaloSubdetectorGeometry &barrelGeom, const CaloSubdetectorTopology &endcapTopo, const CaloSubdetectorGeometry &endcapGeom)
void findRecHitNeighboursCT(reco::PFRecHit &rh, const std::map< unsigned, unsigned > &sortedHits, const CaloSubdetectorTopology &topology)
virtual std::vector< DetId > north(const DetId &id) const =0
int iphi() const
get the cell iphi
CaloTowerDetId id() const
const HcalChannelQuality * theHcalChStatus
void createRecHits(std::vector< reco::PFRecHit > &rechits, std::vector< reco::PFRecHit > &rechitsCleaned, edm::Event &, const edm::EventSetup &)
int hcalHFDigiTimeFlagValue_
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
int hcalHFInTimeWindowFlagValue_
int HcalMaxAllowedHFInTimeWindowSev_
const CaloTowerConstituentsMap * theTowerConstituentsMap
ESHandle< TrackerGeometry > geometry
std::map< unsigned, unsigned >::const_iterator IDH
virtual std::vector< DetId > south(const DetId &id) const =0
const_iterator find(uint32_t rawId) const
double thresh_Barrel_
rechits with E < threshold will not give rise to a PFRecHit
const_iterator end() const
double thresh_HF_
threshold for HF
DetId getNorth(const DetId &id, const CaloSubdetectorTopology &topology)
double maxShortTiming_Cut
uint32_t getValue() const
virtual T south() const
move the navigator south
const Item * getValues(DetId fId) const
virtual std::vector< DetId > east(const DetId &id) const =0
PFRecHitProducerHCAL(const edm::ParameterSet &)
virtual const CornersVec & getCorners() const =0
Returns the corner points of this cell's volume.