CMS 3D CMS Logo
eECALMatrix.cc
Go to the documentation of this file.
1 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
2 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
3 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
4 #include "DataFormats/EcalDetId/interface/EBDetId.h"
5 #include "DataFormats/EcalDetId/interface/EEDetId.h"
6 #include "DataFormats/EcalDetId/interface/ESDetId.h"
7 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
8 #include "FWCore/MessageLogger/interface/MessageLogger.h"
9 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
10 
11 #include <sstream>
12 
13 namespace spr {
14 
15  std::pair<double, bool> energyECAL(const DetId& id,
16  edm::Handle<EcalRecHitCollection>& hitsEC,
17  const EcalSeverityLevelAlgo* sevlv,
18  bool testSpike,
19  double tMin,
20  double tMax,
21  bool debug) {
22  std::vector<EcalRecHitCollection::const_iterator> hits;
23  spr::findHit(hitsEC, id, hits, debug);
24 
25  std::ostringstream st1;
26  if (debug)
27  st1 << "Xtal 0x" << std::hex << id() << std::dec;
28  const EcalRecHitCollection* recHitsEC = (hitsEC.isValid()) ? hitsEC.product() : nullptr;
29  bool flag = (!testSpike) ? true : (sevlv->severityLevel(id, (*recHitsEC)) != EcalSeverityLevel::kWeird);
30  double ener(0);
31  for (const auto& hit : hits) {
32  double en(0), tt(0);
33  if (hit != hitsEC->end()) {
34  en = hit->energy();
35  tt = hit->time();
36  }
37  if (debug)
38  st1 << " " << tt << " " << en;
39  if (tt > tMin && tt < tMax)
40  ener += en;
41  }
42  if (debug) {
43  if (!flag)
44  st1 << " detected to be a spike";
45  edm::LogVerbatim("IsoTrack") << st1.str();
46  }
47  return std::pair<double, bool>(ener, flag);
48  }
49 
50  std::pair<double, bool> energyECAL(const std::vector<DetId>& vdets,
51  edm::Handle<EcalRecHitCollection>& hitsEC,
52  const EcalSeverityLevelAlgo* sevlv,
53  bool noThrCut,
54  bool testSpike,
55  double eThr,
56  double tMin,
57  double tMax,
58  bool debug) {
59  bool flag(true);
60  double energySum(0.0);
61  for (const auto& id : vdets) {
62  if (id != DetId(0)) {
63  std::pair<double, bool> ecalEn = spr::energyECAL(id, hitsEC, sevlv, testSpike, tMin, tMax, debug);
64  if (!ecalEn.second)
65  flag = false;
66  if ((ecalEn.first > eThr) || noThrCut)
67  energySum += ecalEn.first;
68  }
69  }
70  if (debug)
71  edm::LogVerbatim("IsoTrack") << "energyECAL: energySum = " << energySum << " flag = " << flag;
72  return std::pair<double, bool>(energySum, flag);
73  }
74 
75  std::pair<double, bool> eECALmatrix(const DetId& detId,
76  edm::Handle<EcalRecHitCollection>& hitsEB,
77  edm::Handle<EcalRecHitCollection>& hitsEE,
78  const EcalChannelStatus& chStatus,
79  const CaloGeometry* geo,
80  const CaloTopology* caloTopology,
81  const EcalSeverityLevelAlgo* sevlv,
82  int ieta,
83  int iphi,
84  double ebThr,
85  double eeThr,
86  double tMin,
87  double tMax,
88  bool debug) {
89  std::vector<DetId> vdets;
90  spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
91  if (debug)
92  edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
93  << vdets.size();
94  bool flag(true);
95  for (const auto& id : vdets) {
96  if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalBarrel)) {
97  if (sevlv->severityLevel(id, (*hitsEB)) == EcalSeverityLevel::kWeird)
98  flag = false;
99  } else if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalEndcap)) {
100  if (sevlv->severityLevel(id, (*hitsEE)) == EcalSeverityLevel::kWeird)
101  flag = false;
102  }
103  }
104  return std::pair<double, bool>(spr::energyECAL(vdets, hitsEB, hitsEE, ebThr, eeThr, tMin, tMax, debug), flag);
105  }
106 
107  std::pair<double, bool> eECALmatrix(const DetId& detId,
108  edm::Handle<EcalRecHitCollection>& hitsEB,
109  edm::Handle<EcalRecHitCollection>& hitsEE,
110  const EcalChannelStatus& chStatus,
111  const CaloGeometry* geo,
112  const CaloTopology* caloTopology,
113  const EcalSeverityLevelAlgo* sevlv,
114  const EcalTrigTowerConstituentsMap& ttMap,
115  int ieta,
116  int iphi,
117  double ebThr,
118  double eeThr,
119  double tMin,
120  double tMax,
121  bool debug) {
122  std::vector<DetId> vdets;
123  spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
124  if (debug)
125  edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
126  << vdets.size();
127 
128  bool flag(true);
129  double energySum = 0.0;
130  for (const auto& id : vdets) {
131  if ((id != DetId(0)) && (id.det() == DetId::Ecal) &&
132  ((id.subdetId() == EcalBarrel) || (id.subdetId() == EcalEndcap))) {
133  double eTower = spr::energyECALTower(id, hitsEB, hitsEE, ttMap, debug);
134  bool ok = (id.subdetId() == EcalBarrel) ? (eTower > ebThr) : (eTower > eeThr);
135  if (debug && (!ok))
136  edm::LogVerbatim("IsoTrack") << "Crystal 0x" << std::hex << id() << std::dec << " Flag " << ok;
137  if (ok) {
138  std::pair<double, bool> ecalEn = (id.subdetId() == EcalBarrel)
139  ? spr::energyECAL(id, hitsEB, sevlv, true, tMin, tMax, debug)
140  : spr::energyECAL(id, hitsEE, sevlv, false, tMin, tMax, debug);
141  if (!ecalEn.second)
142  flag = false;
143  energySum += ecalEn.first;
144  }
145  }
146  }
147  if (debug)
148  edm::LogVerbatim("IsoTrack") << "energyECAL: energySum = " << energySum << " flag = " << flag;
149  return std::pair<double, bool>(energySum, flag);
150  }
151 
152  std::pair<double, bool> eECALmatrix(const HcalDetId& detId,
153  edm::Handle<EcalRecHitCollection>& hitsEB,
154  edm::Handle<EcalRecHitCollection>& hitsEE,
155  const CaloGeometry* geo,
156  const CaloTowerConstituentsMap* ctmap,
157  const EcalSeverityLevelAlgo* sevlv,
158  double ebThr,
159  double eeThr,
160  double tMin,
161  double tMax,
162  bool debug) {
163  CaloTowerDetId tower = ctmap->towerOf(detId);
164  std::vector<DetId> ids = ctmap->constituentsOf(tower);
165  if (debug) {
166  edm::LogVerbatim("IsoTrack") << "eECALmatrix: " << detId << " belongs to " << tower << " which has " << ids.size()
167  << " constituents";
168  for (unsigned int i = 0; i < ids.size(); ++i) {
169  std::ostringstream st1;
170  st1 << "[" << i << "] " << std::hex << ids[i].rawId() << std::dec;
171  if (ids[i].det() == DetId::Ecal && ids[i].subdetId() == EcalBarrel) {
172  st1 << " " << EBDetId(ids[i]);
173  } else if (ids[i].det() == DetId::Ecal && ids[i].subdetId() == EcalEndcap) {
174  st1 << " " << EEDetId(ids[i]);
175  } else if (ids[i].det() == DetId::Ecal && ids[i].subdetId() == EcalPreshower) {
176  st1 << " " << ESDetId(ids[i]);
177  } else if (ids[i].det() == DetId::Hcal) {
178  st1 << " " << HcalDetId(ids[i]);
179  }
180  edm::LogVerbatim("IsoTrack") << st1.str();
181  }
182  }
183 
184  std::vector<DetId> idEBEE;
185  bool flag(true);
186  for (const auto& id : ids) {
187  if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalBarrel)) {
188  idEBEE.emplace_back(id);
189  if (sevlv->severityLevel(id, (*hitsEB)) == EcalSeverityLevel::kWeird)
190  flag = false;
191  } else if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalEndcap)) {
192  idEBEE.emplace_back(id);
193  if (sevlv->severityLevel(id, (*hitsEE)) == EcalSeverityLevel::kWeird)
194  flag = false;
195  }
196  }
197 
198  if (debug)
199  edm::LogVerbatim("IsoTrack") << "eECALmatrix: with " << idEBEE.size() << " EB+EE hits and "
200  << "spike flag " << flag;
201  double etot = (!idEBEE.empty()) ? spr::energyECAL(idEBEE, hitsEB, hitsEE, ebThr, eeThr, tMin, tMax, debug) : 0;
202  return std::pair<double, bool>(etot, flag);
203  }
204 
205  std::pair<double, bool> eECALmatrix(const DetId& detId,
206  edm::Handle<EcalRecHitCollection>& hitsEB,
207  edm::Handle<EcalRecHitCollection>& hitsEE,
208  const EcalChannelStatus& chStatus,
209  const CaloGeometry* geo,
210  const CaloTopology* caloTopology,
211  const EcalSeverityLevelAlgo* sevlv,
212  const EcalPFRecHitThresholds* eThresholds,
213  int ieta,
214  int iphi,
215  bool debug) {
216  std::vector<DetId> vdets;
217  spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
218  if (debug)
219  edm::LogVerbatim("IsoTrack") << "Inside eECALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " nXtals "
220  << vdets.size();
221  bool flag(true);
222  for (const auto& id : vdets) {
223  if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalBarrel)) {
224  if (sevlv->severityLevel(id, (*hitsEB)) == EcalSeverityLevel::kWeird)
225  flag = false;
226  } else if ((id.det() == DetId::Ecal) && (id.subdetId() == EcalEndcap)) {
227  if (sevlv->severityLevel(id, (*hitsEE)) == EcalSeverityLevel::kWeird)
228  flag = false;
229  }
230  }
231  double energySum = 0.0;
232  for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
233  if (vdets[i1] != DetId(0)) {
234  std::vector<EcalRecHitCollection::const_iterator> hit;
235  if (vdets[i1].subdetId() == EcalBarrel) {
236  spr::findHit(hitsEB, vdets[i1], hit, debug);
237  } else if (vdets[i1].subdetId() == EcalEndcap) {
238  spr::findHit(hitsEE, vdets[i1], hit, debug);
239  }
240  std::ostringstream st1;
241  if (debug)
242  st1 << "Crystal 0x" << std::hex << vdets[i1]() << std::dec;
243  double ener = 0, ethr = static_cast<double>((*eThresholds)[vdets[i1]]);
244  for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
245  double en = 0;
246  if (vdets[i1].subdetId() == EcalBarrel) {
247  if (hit[ihit] != hitsEB->end())
248  en = hit[ihit]->energy();
249  } else if (vdets[i1].subdetId() == EcalEndcap) {
250  if (hit[ihit] != hitsEE->end())
251  en = hit[ihit]->energy();
252  }
253  if (debug)
254  st1 << " " << ihit << " " << en << " Thr " << ethr;
255  ener += en;
256  }
257  if (debug)
258  edm::LogVerbatim("IsoTrack") << st1.str();
259  if (ener > ethr)
260  energySum += ener;
261  }
262  }
263  return std::pair<double, bool>(energySum, flag);
264  }
265 
266 } // namespace spr
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:134
mps_fire.i
i
Definition: mps_fire.py:429
spr::energyECAL
double energyECAL(std::vector< DetId > &vdets, edm::Handle< T > &hitsEB, edm::Handle< T > &hitsEE, double ebThr=-100, double eeThr=-100, double tMin=-500, double tMax=500, bool debug=false)
spr::findHit
std::vector< typename T::const_iterator > findHit(edm::Handle< T > &hits, DetId thisDet, bool debug=false)
spr::matrixECALIds
void matrixECALIds(const DetId &det, int ieta, int iphi, const CaloGeometry *geo, const CaloTopology *caloTopology, std::vector< DetId > &vdets, bool debug=false, bool igNoreTransition=true)
Definition: MatrixECALDetIds.cc:16
CaloTowerConstituentsMap::towerOf
CaloTowerDetId towerOf(const DetId &id) const
Get the tower id for this det id (or null if not known)
Definition: CaloTowerConstituentsMap.cc:27
hfClusterShapes_cfi.hits
hits
Definition: hfClusterShapes_cfi.py:5
edm::Handle::product
T const * product() const
Definition: Handle.h:70
l1ctLayer2EG_cff.id
id
Definition: l1ctLayer2EG_cff.py:85
hcalRecHitTable_cff.energy
energy
Definition: hcalRecHitTable_cff.py:13
EcalSeverityLevelAlgo::severityLevel
EcalSeverityLevel::SeverityLevel severityLevel(const DetId &id) const
Evaluate status from id use channelStatus from DB.
Definition: EcalSeverityLevelAlgo.cc:85
hcalRecHitTable_cff.ieta
ieta
Definition: hcalRecHitTable_cff.py:15
tt
Definition: TTTypes.h:54
hcalRecHitTable_cff.iphi
iphi
Definition: hcalRecHitTable_cff.py:16
hcalRecHitTable_cff.detId
detId
Definition: hcalRecHitTable_cff.py:12
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:267
DetId
Definition: DetId.h:17
debug
#define debug
Definition: HDRShower.cc:19
testProducerWithPsetDescEmpty_cfi.i1
i1
Definition: testProducerWithPsetDescEmpty_cfi.py:45
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:70
CastorDataFrameFilter_impl::energySum
double energySum(const DataFrame &df, int fs, int ls)
Definition: CastorDataFrameFilter.cc:22
spr::energyECALTower
double energyECALTower(const DetId &detId, edm::Handle< T > &hitsEB, edm::Handle< T > &hitsEE, const EcalTrigTowerConstituentsMap &ttMap, bool debug=false)
CaloTowerConstituentsMap::constituentsOf
std::vector< DetId > constituentsOf(const CaloTowerDetId &id) const
Get the constituent detids for this tower id ( not yet implemented )
Definition: CaloTowerConstituentsMap.cc:70
l1tHGCalTowerProducer_cfi.tower
tower
Definition: l1tHGCalTowerProducer_cfi.py:4
spr::eECALmatrix
double eECALmatrix(const DetId &detId, edm::Handle< T > &hitsEB, edm::Handle< T > &hitsEE, const CaloGeometry *geo, const CaloTopology *caloTopology, int ieta, int iphi, double ebThr=-100, double eeThr=-100, double tMin=-500, double tMax=500, bool debug=false)
Generated for CMSSW Reference Manual by doxygen 1.8.14