CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
eECALMatrix.cc
Go to the documentation of this file.
4 #include<iostream>
5 
6 namespace spr{
7 
8  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, double eeThr, double tMin, double tMax, bool debug) {
9 
10  std::vector<DetId> vdets;
11  spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
12 
13  const EcalRecHitCollection * recHitsEB = 0;
14  if (hitsEB.isValid()) recHitsEB = hitsEB.product();
15  bool flag = true;
16  if (debug) {
17  std::cout << "Inside eECALmatrix " << 2*ieta+1 << "X" << 2*iphi+1
18  << " nXtals " << vdets.size() << std::endl;
19  }
20 
21  double energySum = 0.0;
22  for (unsigned int i1=0; i1<vdets.size(); i1++) {
23  if (vdets[i1] != DetId(0)) {
24  bool ok = true;
25  std::vector<EcalRecHitCollection::const_iterator> hit;
26  if (vdets[i1].subdetId()==EcalBarrel) {
27  spr::findHit(hitsEB,vdets[i1],hit,debug);
28 
29  ok = (sevlv->severityLevel(vdets[i1], (*recHitsEB)) != EcalSeverityLevel::kWeird);
30  } else if (vdets[i1].subdetId()==EcalEndcap) {
31  spr::findHit(hitsEE,vdets[i1],hit,debug);
32  }
33  if (debug) std::cout << "Xtal 0x" <<std::hex << vdets[i1]() <<std::dec;
34  double ener=0, ethr=ebThr;
35  if (vdets[i1].subdetId() !=EcalBarrel) ethr = eeThr;
36  for (unsigned int ihit=0; ihit<hit.size(); ihit++) {
37  double en=0, tt=0;
38  if (vdets[i1].subdetId()==EcalBarrel) {
39  if (hit[ihit] != hitsEB->end()) {
40  en = hit[ihit]->energy();
41  tt = hit[ihit]->time();
42  }
43  } else if (vdets[i1].subdetId()==EcalEndcap) {
44  if (hit[ihit] != hitsEE->end()) {
45  en = hit[ihit]->energy();
46  tt = hit[ihit]->time();
47  }
48  }
49  if (debug) std::cout << " " << ihit << " " << en;
50  if (tt > tMin && tt < tMax) ener += en;
51  }
52  if (!ok) {
53  flag = false;
54  if (debug) std::cout << " detected to be a spike";
55  }
56  if (debug) std::cout << "\n";
57  if (ener > ethr) energySum += ener;
58  }
59  }
60  if (debug) std::cout << "energyECAL: energySum = " << energySum << " flag = " << flag << std::endl;
61  return std::pair<double,bool>(energySum,flag);
62  }
63 
64  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, double eeThr, double tMin, double tMax, bool debug) {
65 
66  std::vector<DetId> vdets;
67  spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
68 
69  const EcalRecHitCollection * recHitsEB = 0;
70  if (hitsEB.isValid()) recHitsEB = hitsEB.product();
71  bool flag = true;
72  if (debug) {
73  std::cout << "Inside eECALmatrix " << 2*ieta+1 << "X" << 2*iphi+1
74  << " nXtals " << vdets.size() << std::endl;
75  }
76 
77  double energySum = 0.0;
78  for (unsigned int i1=0; i1<vdets.size(); i1++) {
79  if (vdets[i1] != DetId(0)) {
80  double eTower = spr::energyECALTower(vdets[i1], hitsEB, hitsEE, ttMap, debug);
81  bool ok = true;
82  if (vdets[i1].subdetId()==EcalBarrel) ok = (eTower > ebThr);
83  else if (vdets[i1].subdetId()==EcalEndcap) ok = (eTower > eeThr);
84  if (debug) std::cout << "Crystal 0x" <<std::hex << vdets[i1]()
85  <<std::dec << " Flag " << ok;
86  if (ok) {
87  std::vector<EcalRecHitCollection::const_iterator> hit;
88  if (vdets[i1].subdetId()==EcalBarrel) {
89  spr::findHit(hitsEB,vdets[i1],hit,debug);
90 
91  ok = (sevlv->severityLevel(vdets[i1], (*recHitsEB)) != EcalSeverityLevel::kWeird);
92  } else if (vdets[i1].subdetId()==EcalEndcap) {
93  spr::findHit(hitsEE,vdets[i1],hit,debug);
94  }
95  double ener=0;
96  for (unsigned int ihit=0; ihit<hit.size(); ihit++) {
97  double en=0, tt=0;
98  if (vdets[i1].subdetId()==EcalBarrel) {
99  if (hit[ihit] != hitsEB->end()) {
100  en = hit[ihit]->energy();
101  tt = hit[ihit]->time();
102  }
103  } else if (vdets[i1].subdetId()==EcalEndcap) {
104  if (hit[ihit] != hitsEE->end()) {
105  en = hit[ihit]->energy();
106  tt = hit[ihit]->time();
107  }
108  }
109  if (debug) std::cout << " " << ihit << " E " << en << " T " << tt;
110  if (tt > tMin && tt < tMax) ener += en;
111  }
112  if (!ok) {
113  flag = false;
114  if (debug) std::cout << " detected to be a spike";
115  }
116  energySum += ener;
117  }
118  if (debug) std::cout << "\n";
119  }
120  }
121  if (debug) std::cout << "energyECAL: energySum = " << energySum << " flag = " << flag << std::endl;
122  return std::pair<double,bool>(energySum,flag);
123  }
124 }
EcalSeverityLevel::SeverityLevel severityLevel(const DetId &id) const
Evaluate status from id use channelStatus from DB.
std::vector< typename T::const_iterator > findHit(edm::Handle< T > &hits, DetId thisDet, bool debug=false)
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)
bool isValid() const
Definition: HandleBase.h:75
Definition: DetId.h:18
#define debug
Definition: HDRShower.cc:19
T const * product() const
Definition: Handle.h:81
tuple cout
Definition: gather_cfg.py:121
double energySum(const DataFrame &df, int fs, int ls)
double energyECALTower(const DetId &detId, edm::Handle< T > &hitsEB, edm::Handle< T > &hitsEE, const EcalTrigTowerConstituentsMap &ttMap, bool debug=false)
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)