00001 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
00002 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
00003 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00004 #include<iostream>
00005
00006 namespace spr{
00007
00008 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) {
00009
00010 std::vector<DetId> vdets;
00011 spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
00012
00013 const EcalRecHitCollection * recHitsEB = 0;
00014 if (hitsEB.isValid()) recHitsEB = hitsEB.product();
00015 bool flag = true;
00016 if (debug) {
00017 std::cout << "Inside eECALmatrix " << 2*ieta+1 << "X" << 2*iphi+1
00018 << " nXtals " << vdets.size() << std::endl;
00019 }
00020
00021 double energySum = 0.0;
00022 for (unsigned int i1=0; i1<vdets.size(); i1++) {
00023 if (vdets[i1] != DetId(0)) {
00024 bool ok = true;
00025 std::vector<EcalRecHitCollection::const_iterator> hit;
00026 if (vdets[i1].subdetId()==EcalBarrel) {
00027 spr::findHit(hitsEB,vdets[i1],hit,debug);
00028
00029 ok = (sevlv->severityLevel(vdets[i1], (*recHitsEB)) != EcalSeverityLevel::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 EcalSeverityLevelAlgo* sevlv,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
00091 ok = (sevlv->severityLevel(vdets[i1], (*recHitsEB)) != EcalSeverityLevel::kWeird);
00092 } else if (vdets[i1].subdetId()==EcalEndcap) {
00093 spr::findHit(hitsEE,vdets[i1],hit,debug);
00094 }
00095 double ener=0;
00096 for (unsigned int ihit=0; ihit<hit.size(); ihit++) {
00097 double en=0, tt=0;
00098 if (vdets[i1].subdetId()==EcalBarrel) {
00099 if (hit[ihit] != hitsEB->end()) {
00100 en = hit[ihit]->energy();
00101 tt = hit[ihit]->time();
00102 }
00103 } else if (vdets[i1].subdetId()==EcalEndcap) {
00104 if (hit[ihit] != hitsEE->end()) {
00105 en = hit[ihit]->energy();
00106 tt = hit[ihit]->time();
00107 }
00108 }
00109 if (debug) std::cout << " " << ihit << " E " << en << " T " << tt;
00110 if (tt > tMin && tt < tMax) ener += en;
00111 }
00112 if (!ok) {
00113 flag = false;
00114 if (debug) std::cout << " detected to be a spike";
00115 }
00116 energySum += ener;
00117 }
00118 if (debug) std::cout << "\n";
00119 }
00120 }
00121 if (debug) std::cout << "energyECAL: energySum = " << energySum << " flag = " << flag << std::endl;
00122 return std::pair<double,bool>(energySum,flag);
00123 }
00124 }
1.7.3