Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Run2-alca94 Add a new method to compute energy sum of ECAL towers preceding a HCAl cell… #20427

Merged
merged 2 commits into from Sep 12, 2017
Merged

Run2-alca94 Add a new method to compute energy sum of ECAL towers preceding a HCAl cell… #20427

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a3d5943
Add a new method to compute energy sum of ECAL towers preceding a HCA…
Sep 7, 2017
f9d6053
Simplify eECALMatrix a bit
Sep 11, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
4 changes: 4 additions & 0 deletions Calibration/IsolatedParticles/interface/eECALMatrix.h
View file
Open in desktop
Expand Up @@ -27,6 +27,7 @@ Created: August 2009
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"

Expand All @@ -35,6 +36,7 @@ Created: August 2009
#include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloTopology/interface/CaloTopology.h"
#include "Geometry/CaloTopology/interface/CaloTowerConstituentsMap.h"
#include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"

class EcalSeverityLevelAlgo;
Expand All @@ -54,6 +56,8 @@ namespace spr{
std::pair <double,bool> eECALmatrix(const DetId& detId, edm::Handle<EcalRecHitCollection>& hitsEB, edm::Handle<EcalRecHitCollection>& hitsEE, const EcalChannelStatus& chStatus, const CaloGeometry* geo, const CaloTopology* caloTopology, const EcalSeverityLevelAlgo* sevlv,int ieta, int iphi, double ebThr=-100, double eeThr=-100, double tMin=-500, double tMax=500, bool debug=false);

std::pair <double,bool> eECALmatrix(const DetId& detId, edm::Handle<EcalRecHitCollection>& hitsEB, edm::Handle<EcalRecHitCollection>& hitsEE, const EcalChannelStatus& chStatus, const CaloGeometry* geo, const CaloTopology* caloTopology, const EcalSeverityLevelAlgo* sevlv,const EcalTrigTowerConstituentsMap& ttMap, int ieta, int iphi, double ebThr=-100, double eeThr=-100, double tMin=-500, double tMax=500, bool debug=false);

std::pair<double,bool> eECALmatrix(const HcalDetId& detId, edm::Handle<EcalRecHitCollection>& hitsEB, edm::Handle<EcalRecHitCollection>& hitsEE, const CaloGeometry* geo, const CaloTowerConstituentsMap* ctmap, const EcalSeverityLevelAlgo* sevlv, double ebThr=-100, double eeThr=-100, double tMin=-500, double tMax=500, bool debug=false);

// returns vector of hits in NxN matrix
template <typename T>
Expand Down
238 changes: 130 additions & 108 deletions Calibration/IsolatedParticles/src/eECALMatrix.cc
View file
Open in desktop
@@ -1,12 +1,73 @@
#include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
#include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
#include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/EcalDetId/interface/ESDetId.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
#include<iostream>

//#define EDM_ML_DEBUG

namespace spr{

std::pair<double,bool> energyECAL(const DetId& id,
edm::Handle<EcalRecHitCollection>& hitsEC,
const EcalSeverityLevelAlgo* sevlv,
bool testSpike, double tMin, double tMax,
bool debug) {
std::vector<EcalRecHitCollection::const_iterator> hits;
spr::findHit(hitsEC,id,hits,debug);
#ifdef EDM_ML_DEBUG
if (debug) std::cout << "Xtal 0x" << std::hex << id() <<std::dec;
#endif
const EcalRecHitCollection* recHitsEC = (hitsEC.isValid()) ? hitsEC.product() : nullptr;
bool flag = (!testSpike) ? true :
(sevlv->severityLevel(id,(*recHitsEC)) != EcalSeverityLevel::kWeird);
double ener(0);
for (const auto& hit : hits) {
double en(0), tt(0);
if (hit != hitsEC->end()) {
en = hit->energy();
tt = hit->time();
}
#ifdef EDM_ML_DEBUG
if (debug) std::cout << " " << tt << " " << en;
#endif
if (tt > tMin && tt < tMax) ener += en;
}
#ifdef EDM_ML_DEBUG
if (!flag && debug) std::cout << " detected to be a spike";
if (debug) std::cout << std::endl;
#endif
return std::pair<double,bool>(ener,flag);
}

std::pair<double,bool> energyECAL(const std::vector<DetId>& vdets,
edm::Handle<EcalRecHitCollection>& hitsEC,
const EcalSeverityLevelAlgo* sevlv,
bool noThrCut, bool testSpike, double eThr,
double tMin, double tMax, bool debug) {

bool flag(true);
double energySum(0.0);
for (const auto& id : vdets) {
if (id != DetId(0)) {
std::pair<double,bool> ecalEn = spr::energyECAL(id,hitsEC,sevlv,
testSpike,tMin,tMax,
debug);
if (!ecalEn.second) flag = false;
if ((ecalEn.first>eThr) || noThrCut) energySum += ecalEn.first;
}
}
#ifdef EDM_ML_DEBUG
if (debug) std::cout << "energyECAL: energySum = " << energySum
<< " flag = " << flag << std::endl;
#endif
return std::pair<double,bool>(energySum,flag);
}

std::pair<double,bool> eECALmatrix(const DetId& detId,
edm::Handle<EcalRecHitCollection>& hitsEB,
edm::Handle<EcalRecHitCollection>& hitsEE,
Expand All @@ -17,72 +78,25 @@ namespace spr{
int ieta, int iphi, double ebThr,
double eeThr, double tMin, double tMax,
bool debug) {

std::vector<DetId> vdets;
spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);

const EcalRecHitCollection * recHitsEB = 0;
if (hitsEB.isValid()) recHitsEB = hitsEB.product();
bool flag = true;
#ifdef EDM_ML_DEBUG
if (debug) {
std::cout << "Inside eECALmatrix " << 2*ieta+1 << "X" << 2*iphi+1
<< " nXtals " << vdets.size() << std::endl;
}
#endif
double energySum = 0.0;
for (unsigned int i1=0; i1<vdets.size(); i1++) {
if (vdets[i1] != DetId(0)) {
bool ok = true;
std::vector<EcalRecHitCollection::const_iterator> hit;
if (vdets[i1].subdetId()==EcalBarrel) {
spr::findHit(hitsEB,vdets[i1],hit,debug);

ok = (sevlv->severityLevel(vdets[i1], (*recHitsEB)) != EcalSeverityLevel::kWeird);
} else if (vdets[i1].subdetId()==EcalEndcap) {
spr::findHit(hitsEE,vdets[i1],hit,debug);
}
#ifdef EDM_ML_DEBUG
if (debug) std::cout << "Xtal 0x" <<std::hex << vdets[i1]() <<std::dec;
#endif
double ener=0, ethr=ebThr;
if (vdets[i1].subdetId() !=EcalBarrel) ethr = eeThr;
for (unsigned int ihit=0; ihit<hit.size(); ihit++) {
double en=0, tt=0;
if (vdets[i1].subdetId()==EcalBarrel) {
if (hit[ihit] != hitsEB->end()) {
en = hit[ihit]->energy();
tt = hit[ihit]->time();
}
} else if (vdets[i1].subdetId()==EcalEndcap) {
if (hit[ihit] != hitsEE->end()) {
en = hit[ihit]->energy();
tt = hit[ihit]->time();
}
}
#ifdef EDM_ML_DEBUG
if (debug) std::cout << " " << ihit << " " << en;
#endif
if (tt > tMin && tt < tMax) ener += en;
}
if (!ok) {
flag = false;
#ifdef EDM_ML_DEBUG
if (debug) std::cout << " detected to be a spike";
#endif
}
#ifdef EDM_ML_DEBUG
if (debug) std::cout << std::endl;
#endif
if (ener > ethr) energySum += ener;
}
}
#ifdef EDM_ML_DEBUG
if (debug) std::cout << "energyECAL: energySum = " << energySum << " flag = " << flag << std::endl;
#endif
return std::pair<double,bool>(energySum,flag);

if (detId.det() == DetId::Ecal && detId.subdetId() == EcalBarrel) {
return spr::energyECAL(vdets,hitsEB,sevlv,false,true,ebThr,tMin,tMax,debug);
} else if (detId.det() == DetId::Ecal && detId.subdetId() == EcalEndcap) {
return spr::energyECAL(vdets,hitsEE,sevlv,false,false,eeThr,tMin,tMax,debug);
} else {
return std::pair<double,bool>(0,true);
}
}

std::pair<double,bool> eECALmatrix(const DetId& detId,
edm::Handle<EcalRecHitCollection>& hitsEB,
edm::Handle<EcalRecHitCollection>& hitsEE,
Expand All @@ -94,74 +108,82 @@ namespace spr{
int ieta, int iphi, double ebThr,
double eeThr, double tMin, double tMax,
bool debug) {

std::vector<DetId> vdets;
spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);

const EcalRecHitCollection * recHitsEB = 0;
if (hitsEB.isValid()) recHitsEB = hitsEB.product();
bool flag = true;
#ifdef EDM_ML_DEBUG
if (debug) {
std::cout << "Inside eECALmatrix " << 2*ieta+1 << "X" << 2*iphi+1
<< " nXtals " << vdets.size() << std::endl;
}
}
#endif

bool flag(true);
double energySum = 0.0;
for (unsigned int i1=0; i1<vdets.size(); i1++) {
if (vdets[i1] != DetId(0)) {
double eTower = spr::energyECALTower(vdets[i1], hitsEB, hitsEE, ttMap, debug);
bool ok = true;
if (vdets[i1].subdetId()==EcalBarrel) ok = (eTower > ebThr);
else if (vdets[i1].subdetId()==EcalEndcap) ok = (eTower > eeThr);
for (const auto & id : vdets) {
if ((id != DetId(0)) && (id.det() == DetId::Ecal) &&
((id.subdetId()==EcalBarrel) || (id.subdetId()==EcalEndcap))) {
double eTower = spr::energyECALTower(id, hitsEB, hitsEE, ttMap, debug);
bool ok = (id.subdetId()==EcalBarrel) ? (eTower > ebThr) : (eTower > eeThr);
#ifdef EDM_ML_DEBUG
if (debug) std::cout << "Crystal 0x" <<std::hex << vdets[i1]()
<<std::dec << " Flag " << ok;
if (debug && (!ok)) std::cout << "Crystal 0x" << std::hex << id()
<< std::dec << " Flag " << ok <<std::endl;
#endif
if (ok) {
std::vector<EcalRecHitCollection::const_iterator> hit;
if (vdets[i1].subdetId()==EcalBarrel) {
spr::findHit(hitsEB,vdets[i1],hit,debug);

ok = (sevlv->severityLevel(vdets[i1], (*recHitsEB)) != EcalSeverityLevel::kWeird);
} else if (vdets[i1].subdetId()==EcalEndcap) {
spr::findHit(hitsEE,vdets[i1],hit,debug);
}
double ener=0;
for (unsigned int ihit=0; ihit<hit.size(); ihit++) {
double en=0, tt=0;
if (vdets[i1].subdetId()==EcalBarrel) {
if (hit[ihit] != hitsEB->end()) {
en = hit[ihit]->energy();
tt = hit[ihit]->time();
}
} else if (vdets[i1].subdetId()==EcalEndcap) {
if (hit[ihit] != hitsEE->end()) {
en = hit[ihit]->energy();
tt = hit[ihit]->time();
}
}
std::pair<double,bool> ecalEn = (id.subdetId()==EcalBarrel) ?
spr::energyECAL(id,hitsEB,sevlv,true,tMin,tMax,debug) :
spr::energyECAL(id,hitsEE,sevlv,false,tMin,tMax,debug);
if (!ecalEn.second) flag = false;
energySum += ecalEn.first;
}
}
}
#ifdef EDM_ML_DEBUG
if (debug) std::cout << " " << ihit << " E " << en << " T " << tt;
if (debug) std::cout << "energyECAL: energySum = " << energySum
<< " flag = " << flag << std::endl;
#endif
if (tt > tMin && tt < tMax) ener += en;
}
if (!ok) {
flag = false;
return std::pair<double,bool>(energySum,flag);
}

std::pair<double,bool> eECALmatrix(const HcalDetId& detId,
edm::Handle<EcalRecHitCollection>& hitsEB,
edm::Handle<EcalRecHitCollection>& hitsEE,
const CaloGeometry* geo,
const CaloTowerConstituentsMap* ctmap,
const EcalSeverityLevelAlgo* sevlv,
double ebThr, double eeThr, double tMin,
double tMax, bool debug) {

CaloTowerDetId tower = ctmap->towerOf(detId);
std::vector<DetId> ids = ctmap->constituentsOf(tower);
#ifdef EDM_ML_DEBUG
if (debug) std::cout << " detected to be a spike";
#endif
}
energySum += ener;
if (debug) {
std::cout << "eECALmatrix: " << detId << " belongs to " << tower
<< " which has " << ids.size() << " constituents" << std::endl;
for (unsigned int i=0; i<ids.size(); ++i) {
std::cout << "[" << i << "] " << std::hex << ids[i].rawId() <<std::dec;
if (ids[i].det()==DetId::Ecal && ids[i].subdetId()==EcalBarrel){
std::cout << " " << EBDetId(ids[i]) << std::endl;
} else if (ids[i].det()==DetId::Ecal && ids[i].subdetId()==EcalEndcap){
std::cout << " " << EEDetId(ids[i]) << std::endl;
} else if (ids[i].det()==DetId::Ecal && ids[i].subdetId()==EcalPreshower) {
std::cout << " " << ESDetId(ids[i]) << std::endl;
} else if (ids[i].det()==DetId::Hcal) {
std::cout << " " << HcalDetId(ids[i]) << std::endl;
} else {
std::cout << std::endl;
}
#ifdef EDM_ML_DEBUG
if (debug) std::cout << std::endl;
#endif
}
}
#ifdef EDM_ML_DEBUG
if (debug) std::cout << "energyECAL: energySum = " << energySum << " flag = " << flag << std::endl;
#endif
return std::pair<double,bool>(energySum,flag);

if (detId.det() == DetId::Ecal && detId.subdetId() == EcalBarrel) {
return spr::energyECAL(ids,hitsEB,sevlv,false,true,ebThr,tMin,tMax,debug);
} else if (detId.det() == DetId::Ecal && detId.subdetId() == EcalEndcap) {
return spr::energyECAL(ids,hitsEE,sevlv,false,false,eeThr,tMin,tMax,debug);
} else {
return std::pair<double,bool>(0,true);
}
}
}