00001 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
00002 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
00003 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00004
00005 #include<iostream>
00006
00007 namespace spr{
00008
00009 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, int ieta, int iphi, double ebThr, double eeThr, double tMin, double tMax, bool debug) {
00010
00011 std::vector<DetId> vdets;
00012 spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
00013
00014 const EcalRecHitCollection * recHitsEB = 0;
00015 if (hitsEB.isValid()) recHitsEB = hitsEB.product();
00016 bool flag = true;
00017 if (debug) {
00018 std::cout << "Inside eECALmatrix " << 2*ieta+1 << "X" << 2*iphi+1
00019 << " nXtals " << vdets.size() << std::endl;
00020 }
00021
00022 double energySum = 0.0;
00023 for (unsigned int i1=0; i1<vdets.size(); i1++) {
00024 if (vdets[i1] != DetId(0)) {
00025 bool ok = true;
00026 std::vector<EcalRecHitCollection::const_iterator> hit;
00027 if (vdets[i1].subdetId()==EcalBarrel) {
00028 spr::findHit(hitsEB,vdets[i1],hit,debug);
00029 ok = (EcalSeverityLevelAlgo::severityLevel(vdets[i1], (*recHitsEB), chStatus) != EcalSeverityLevelAlgo::kWeird);
00030 } else if (vdets[i1].subdetId()==EcalEndcap) {
00031 spr::findHit(hitsEE,vdets[i1],hit,debug);
00032 }
00033 if (debug) std::cout << "Xtal 0x" <<std::hex << vdets[i1]() <<std::dec;
00034 double ener=0, ethr=ebThr;
00035 if (vdets[i1].subdetId() !=EcalBarrel) ethr = eeThr;
00036 for (unsigned int ihit=0; ihit<hit.size(); ihit++) {
00037 double en=0, tt=0;
00038 if (vdets[i1].subdetId()==EcalBarrel) {
00039 if (hit[ihit] != hitsEB->end()) {
00040 en = hit[ihit]->energy();
00041 tt = hit[ihit]->time();
00042 }
00043 } else if (vdets[i1].subdetId()==EcalEndcap) {
00044 if (hit[ihit] != hitsEE->end()) {
00045 en = hit[ihit]->energy();
00046 tt = hit[ihit]->time();
00047 }
00048 }
00049 if (debug) std::cout << " " << ihit << " " << en;
00050 if (tt > tMin && tt < tMax) ener += en;
00051 }
00052 if (!ok) {
00053 flag = false;
00054 if (debug) std::cout << " detected to be a spike";
00055 }
00056 if (debug) std::cout << "\n";
00057 if (ener > ethr) energySum += ener;
00058 }
00059 }
00060 if (debug) std::cout << "energyECAL: energySum = " << energySum << " flag = " << flag << std::endl;
00061 return std::pair<double,bool>(energySum,flag);
00062 }
00063
00064 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 EcalTrigTowerConstituentsMap& ttMap, int ieta, int iphi, double ebThr, double eeThr, double tMin, double tMax, bool debug) {
00065
00066 std::vector<DetId> vdets;
00067 spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
00068
00069 const EcalRecHitCollection * recHitsEB = 0;
00070 if (hitsEB.isValid()) recHitsEB = hitsEB.product();
00071 bool flag = true;
00072 if (debug) {
00073 std::cout << "Inside eECALmatrix " << 2*ieta+1 << "X" << 2*iphi+1
00074 << " nXtals " << vdets.size() << std::endl;
00075 }
00076
00077 double energySum = 0.0;
00078 for (unsigned int i1=0; i1<vdets.size(); i1++) {
00079 if (vdets[i1] != DetId(0)) {
00080 double eTower = spr::energyECALTower(vdets[i1], hitsEB, hitsEE, ttMap, debug);
00081 bool ok = true;
00082 if (vdets[i1].subdetId()==EcalBarrel) ok = (eTower > ebThr);
00083 else if (vdets[i1].subdetId()==EcalEndcap) ok = (eTower > eeThr);
00084 if (debug) std::cout << "Crystal 0x" <<std::hex << vdets[i1]()
00085 <<std::dec << " Flag " << ok;
00086 if (ok) {
00087 std::vector<EcalRecHitCollection::const_iterator> hit;
00088 if (vdets[i1].subdetId()==EcalBarrel) {
00089 spr::findHit(hitsEB,vdets[i1],hit,debug);
00090 ok = (EcalSeverityLevelAlgo::severityLevel(vdets[i1], (*recHitsEB), chStatus) != EcalSeverityLevelAlgo::kWeird);
00091 } else if (vdets[i1].subdetId()==EcalEndcap) {
00092 spr::findHit(hitsEE,vdets[i1],hit,debug);
00093 }
00094 double ener=0;
00095 for (unsigned int ihit=0; ihit<hit.size(); ihit++) {
00096 double en=0, tt=0;
00097 if (vdets[i1].subdetId()==EcalBarrel) {
00098 if (hit[ihit] != hitsEB->end()) {
00099 en = hit[ihit]->energy();
00100 tt = hit[ihit]->time();
00101 }
00102 } else if (vdets[i1].subdetId()==EcalEndcap) {
00103 if (hit[ihit] != hitsEE->end()) {
00104 en = hit[ihit]->energy();
00105 tt = hit[ihit]->time();
00106 }
00107 }
00108 if (debug) std::cout << " " << ihit << " E " << en << " T " << tt;
00109 if (tt > tMin && tt < tMax) ener += en;
00110 }
00111 if (!ok) {
00112 flag = false;
00113 if (debug) std::cout << " detected to be a spike";
00114 }
00115 energySum += ener;
00116 }
00117 if (debug) std::cout << "\n";
00118 }
00119 }
00120 if (debug) std::cout << "energyECAL: energySum = " << energySum << " flag = " << flag << std::endl;
00121 return std::pair<double,bool>(energySum,flag);
00122 }
00123 }
1.7.3