![]() |
![]() |
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 }