32 return Calculate(TheCaloGeometry, TheHBHERecHits, TheCaloTowers,TheEBRecHits,TheEERecHits,TheSetup);
40 const int iEtaOverlap = 22;
47 float MinTimeHits[nPhiMax];
49 float MaxTimeHits[nPhiMax];
50 for(
unsigned int i = 0 ;
i < nPhiMax ;
i++ )
58 for (
const auto &
hit : (*TheHBHERecHits)) {
60 switch (
id.subdet() ) {
73 if(iPhi < nPhiMax &&
std::abs(iEta) <= iEtaOverlap) {
74 SumE[iPhi]+=
hit.energy();
78 MinTimeHits[iPhi] = time < MinTimeHits[iPhi] ? time : MinTimeHits[iPhi];
79 MaxTimeHits[iPhi] = time > MaxTimeHits[iPhi] ? time : MaxTimeHits[iPhi];
83 for (
int iPhi = 1 ; iPhi < nPhiMax ; iPhi++ ) {
86 PhiWedge wedge(SumE[iPhi], iPhi, NumHits[iPhi], MinTimeHits[iPhi], MaxTimeHits[iPhi]);
89 std::vector<const HBHERecHit*> Hits;
90 for (
const auto &
hit : (*TheHBHERecHits)) {
92 if(
id.iphi() != iPhi )
continue;
93 if(
std::abs(
id.ieta() ) > iEtaOverlap )
continue;
94 switch (
id.subdet() ) {
104 Hits.push_back(&(
hit));
107 std::sort( Hits.begin() , Hits.end() ,
CompareTime);
108 float MinusToPlus = 0.;
109 float PlusToMinus = 0.;
110 for(
unsigned int i = 0 ;
i < Hits.size() ;
i++ ) {
112 int ieta_i = id_i.
ieta();
113 for(
unsigned int j = (
i+1) ; j < Hits.size() ; j++ ) {
115 int ieta_j = id_j.
ieta();
116 if( ieta_i > ieta_j ) PlusToMinus +=
std::abs(ieta_i - ieta_j ) ;
117 else MinusToPlus +=
std::abs(ieta_i - ieta_j);
120 float PlusZOriginConfidence = (PlusToMinus + MinusToPlus )? PlusToMinus / ( PlusToMinus + MinusToPlus ) : -1. ;
131 std::map<int, float> iPhiHadEtMap;
132 std::vector<const CaloTower*> sortedCaloTowers;
133 for (
const auto & tower : (*TheCaloTowers)) {
134 if(
std::abs(tower.ieta()) > maxAbsIEta)
continue;
136 int iPhi = tower.iphi();
137 if(!iPhiHadEtMap.count(iPhi)) iPhiHadEtMap[iPhi] = 0.0;
138 iPhiHadEtMap[iPhi] += tower.hadEt();
140 if(tower.numProblematicHcalCells() > 0) sortedCaloTowers.push_back(&(tower));
147 std::sort(sortedCaloTowers.begin(), sortedCaloTowers.end(),
CompareTowers);
152 int prevIEta = -99, prevIPhi = -99;
153 float prevHadEt = 0.;
155 std::pair<uint8_t, CaloTowerDetId> prevPair, towerPair;
156 bool wasContiguous =
true;
159 for(
unsigned int i = 0;
i < sortedCaloTowers.size();
i++) {
164 bool newIPhi = tower->
iphi() != prevIPhi;
165 bool isContiguous = tower->
ieta() == 1 ? tower->
ieta() - 2 == prevIEta : tower->
ieta() - 1 == prevIEta;
167 isContiguous = isContiguous || (tower->
ieta() == -maxAbsIEta);
168 if(newIPhi) isContiguous =
false;
170 if(!wasContiguous && isContiguous) {
171 strip.cellTowerIds.push_back(prevPair);
172 strip.cellTowerIds.push_back(towerPair);
173 strip.hadEt += prevHadEt + tower->
hadEt();
174 strip.emEt += prevEmEt + tower->
emEt();
177 if(wasContiguous && isContiguous) {
178 strip.cellTowerIds.push_back(towerPair);
179 strip.hadEt += tower->
hadEt();
180 strip.emEt += tower->
emEt();
183 if((wasContiguous && !isContiguous) ||
i == sortedCaloTowers.size()-1) {
185 if(strip.cellTowerIds.size() > 3) {
187 int iPhi = strip.cellTowerIds.at(0).second.iphi();
188 int iPhiLower = (iPhi == 1) ? 72 : iPhi - 1;
189 int iPhiUpper = (iPhi == 72) ? 1 : iPhi + 1;
191 float energyRatio = 0.0;
192 if(iPhiHadEtMap.count(iPhiLower)) energyRatio += iPhiHadEtMap[iPhiLower];
193 if(iPhiHadEtMap.count(iPhiUpper)) energyRatio += iPhiHadEtMap[iPhiUpper];
194 iPhiHadEtMap[iPhi] =
std::max(iPhiHadEtMap[iPhi], 0.001
F);
196 energyRatio /= iPhiHadEtMap[iPhi];
197 strip.energyRatio = energyRatio;
205 wasContiguous = isContiguous;
206 prevPair = towerPair;
207 prevEmEt = tower->
emEt();
208 prevIPhi = tower->
iphi();
209 prevIEta = tower->
ieta();
210 prevHadEt = tower->
hadEt();
228 std::vector<HaloClusterCandidateHCAL> haloclustercands_HB;
231 std::vector<HaloClusterCandidateHCAL> haloclustercands_HE;
238 return TheHcalHaloData;
251 std::vector<HaloClusterCandidateHCAL> TheHaloClusterCandsHB;
255 for(
size_t ihit = 0; ihit<hbherechitcoll->
size(); ++ ihit){
258 const HBHERecHit & rechit = (*hbherechitcoll)[ ihit ];
261 double rhet = rechit.
energy()*
sqrt(rhpos.perp2()/rhpos.mag2());
262 if(rhet<et_thresh_seedrh)
continue;
263 if(
std::abs(rhpos.z())>zseparation_HBHE)
continue;
264 double eta = rhpos.eta();
265 double phi = rhpos.phi();
269 int nbtowerssameeta(0);
270 double timediscriminatorITBH(0),timediscriminatorOTBH(0);
271 double etstrip_phiseedplus1(0), etstrip_phiseedminus1(0);
275 for(
size_t jhit = 0; jhit<hbherechitcoll->
size(); ++ jhit){
276 const HBHERecHit & rechitj = (*hbherechitcoll)[ jhit ];
279 double rhetj = rechitj.
energy()*
sqrt(rhposj.perp2()/rhposj.mag2());
280 if(rhetj<2)
continue;
281 if(
std::abs(rhposj.z())>zseparation_HBHE)
continue;
282 double etaj = rhposj.eta();
283 double phij = rhposj.phi();
284 double deta = eta - etaj;
290 if(
std::abs(dphi)<0.05) nbtowerssameeta++;
291 if(dphi>0.05) etstrip_phiseedplus1+=rhetj;
292 if(dphi<-0.05) etstrip_phiseedminus1+=rhetj;
305 double rhtj = rechitj.
time();
306 timediscriminatorITBH+= std::log10( rhetj )* ( rhtj +0.5*(
sqrt(240.*240.+rhposj.z()*rhposj.z()) -
std::abs(rhposj.z()))/c_cm_per_ns);
307 if(
std::abs(rhposj.z())<300) timediscriminatorOTBH+= std::log10( rhetj )* ( rhtj -0.5*(25-(
sqrt(330.*330.+rhposj.z()*rhposj.z()) +
std::abs(rhposj.z()))/c_cm_per_ns) );
312 if(etstrip_phiseedplus1/etcluster>0.2&& etstrip_phiseedminus1/etcluster>0.2)
continue;
316 for(
size_t jhit = 0; jhit<ecalrechitcoll->
size(); ++ jhit){
317 const EcalRecHit & rechitj = (*ecalrechitcoll)[ jhit ];
319 double rhetj = rechitj.
energy()*
sqrt(rhposj.perp2()/rhposj.mag2());
320 if(rhetj<2)
continue;
321 double etaj = rhposj.eta();
322 double phij = rhposj.phi();
323 if(
std::abs(eta-etaj)>0.2)
continue;
325 eoh+=rhetj/etcluster;
328 if(eoh>0.1)
continue;
352 TheHaloClusterCandsHB.push_back(clustercand);
355 return TheHaloClusterCandsHB;
361 std::vector<HaloClusterCandidateHCAL> TheHaloClusterCandsHE;
365 for(
size_t ihit = 0; ihit<hbherechitcoll->
size(); ++ ihit){
368 const HBHERecHit & rechit = (*hbherechitcoll)[ ihit ];
371 double rhet = rechit.
energy()*
sqrt(rhpos.perp2()/rhpos.mag2());
372 if(rhet<et_thresh_seedrh)
continue;
373 if(
std::abs(rhpos.z())<zseparation_HBHE)
continue;
374 double eta = rhpos.eta();
375 double phi = rhpos.phi();
376 double rhr =
sqrt(rhpos.perp2());
378 double etcluster(0),hdepth1(0);
380 double etstrip_phiseedplus1(0), etstrip_phiseedminus1(0);
384 for(
size_t jhit = 0; jhit<hbherechitcoll->
size(); ++ jhit){
385 const HBHERecHit & rechitj = (*hbherechitcoll)[ jhit ];
388 double rhetj = rechitj.
energy()*
sqrt(rhposj.perp2()/rhposj.mag2());
389 if(rhetj<2)
continue;
390 if(
std::abs(rhposj.z())<zseparation_HBHE)
continue;
391 if(rhpos.z()*rhposj.z()<0)
continue;
392 double phij = rhposj.phi();
395 double rhrj =
sqrt(rhposj.perp2());
396 if(
std::abs( rhr-rhrj )>50)
continue;
398 if(dphi>0.05) etstrip_phiseedplus1+=rhetj;
399 if(dphi<-0.05) etstrip_phiseedminus1+=rhetj;
402 if(
std::abs( rhposj.z())<405 )hdepth1+=rhetj;
408 if(etstrip_phiseedplus1/etcluster>0.1&& etstrip_phiseedminus1/etcluster>0.1)
continue;
412 for(
size_t jhit = 0; jhit<ecalrechitcoll->
size(); ++ jhit){
413 const EcalRecHit & rechitj = (*ecalrechitcoll)[ jhit ];
415 double rhetj = rechitj.
energy()*
sqrt(rhposj.perp2()/rhposj.mag2());
416 if(rhetj<2)
continue;
417 if(rhpos.z()*rhposj.z()<0)
continue;
418 double etaj = rhposj.eta();
419 double phij = rhposj.phi();
423 eoh+=rhetj/etcluster;
426 if(eoh>0.1)
continue;
450 TheHaloClusterCandsHE.push_back(clustercand);
453 return TheHaloClusterCandsHE;
491 if(hcand.
getSeedR()>170)
return false;
498 if(ishlt && hcand.
getSeedEt()<50)
return false;
void setIsHaloFromPattern(bool x)
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
double getEtStripPhiSeedPlus1() const
float HERecHitEnergyThreshold
void setHaloClusterCandidatesHB(const std::vector< HaloClusterCandidateHCAL > &x)
HcalDetId id() const
get the id
Global3DPoint GlobalPoint
void setH1overH123(double x)
double getH1overH123() const
std::vector< reco::HaloClusterCandidateHCAL > GetHaloClusterCandidateHE(edm::Handle< EcalRecHitCollection > &eerechitcoll, edm::Handle< HBHERecHitCollection > &hbherechitcoll, float et_thresh_seedrh)
const CaloGeometry * geo_
int getNbTowersInEta() const
double getTimeDiscriminatorOTBH() const
bool HBClusterShapeandTimeStudy(reco::HaloClusterCandidateHCAL hcand, bool ishlt)
void setSeedEta(double x)
void setEtStripPhiSeedPlus1(double x)
double getEtStripPhiSeedMinus1() const
int ieta() const
get the cell ieta
void setHaloClusterCandidatesHE(const std::vector< HaloClusterCandidateHCAL > &x)
Abs< T >::type abs(const T &t)
void SetPlusZOriginConfidence(float x)
void setNbTowersInEta(double x)
math::XYZPoint Point
point in the space
const GlobalPoint & getPosition(const DetId &id) const
Get the position of a given detector id.
bool CompareTime(const HBHERecHit *x, const HBHERecHit *y)
GlobalPoint getPosition(const DetId &id) const
double getTimeDiscriminatorITBH() const
void setTimeDiscriminator(double x)
void setSeedTime(double x)
reco::HcalHaloData Calculate(const CaloGeometry &TheCaloGeometry, edm::Handle< HBHERecHitCollection > &TheHBHERecHits, edm::Handle< CaloTowerCollection > &TheCaloTowers, edm::Handle< EBRecHitCollection > &TheEBRecHits, edm::Handle< EERecHitCollection > &TheEERecHits, const edm::EventSetup &TheSetup)
CaloTowerDetId id() const
unsigned int numProblematicHcalCells() const
DetId id() const
get the id
void setEtStripPhiSeedMinus1(double x)
XYZPointD XYZPoint
point in space with cartesian internal representation
bool CompareTowers(const CaloTower *x, const CaloTower *y)
float HBRecHitEnergyThreshold
const HcalGeometry * hgeo_
void setTimeDiscriminatorOTBH(double x)
math::XYZPoint getPosition(const DetId &id, reco::Vertex::Point vtx)
const std::vector< HaloTowerStrip > & getProblematicStrips() const
void setBeamHaloRecHitsCandidates(edm::RefVector< HBHERecHitCollection > x)
void setSeedPhi(double x)
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
const std::vector< PhiWedge > & GetPhiWedges() const
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
void setIsHaloFromPattern_HLT(bool x)
void setClusterSize(int x)
T const * product() const
bool HEClusterShapeandTimeStudy(reco::HaloClusterCandidateHCAL hcand, bool ishlt)
void setTimeDiscriminatorITBH(double x)
void setClusterEt(double x)
std::vector< reco::HaloClusterCandidateHCAL > GetHaloClusterCandidateHB(edm::Handle< EcalRecHitCollection > &ebrechitcoll, edm::Handle< HBHERecHitCollection > &hbherechitcoll, float et_thresh_seedrh)